Variational autoencoder with weighted samples for high-dimensional non-parametric adaptive importance sampling
Authors: Julien Demange-Chryst
Francois Bachoc
Jérôme Morio
Timothé Krauth
Venue: Accepted by TMLR ()

Abstract: Adaptive importance sampling is a well-known family of algorithms for density approximation, generation and Monte Carlo integration including rare event estimation. The main common denominator of this family of algorithms is to perform density estimation with weighted samples at each iteration. However, the classical existing methods to do so, such as kernel smoothing or approximation by a Gaussian distribution, suffer from the curse of dimensionality and/or a lack of flexibility. Both are limitations in high dimension and when we do not have any prior knowledge on the form of the target distribution, such as its number of modes. Variational autoencoders are probabilistic tools able to represent with fidelity high-dimensional data in a lower dimensional space. They constitute a parametric family of distributions robust faced to the dimension and since they are based on deep neural networks, they are flexible enough to be considered as non-parametric models. In this paper, we propose to use a variational autoencoder as the auxiliary importance sampling distribution by extending the existing framework to weighted samples. We integrate the proposed procedure in existing adaptive importance sampling algorithms and we illustrate its practical interest on diverse examples.

Published in Transactions on Machine Learning Research (01/2024)
Variational autoencoder with weighted samples for high- dimensional non-parametric adaptive importance sampling
Julien Demange- Chryst.
juien..emmngg-chrgtffoonne..
ONERA/DTIS, Université de Toulouse, F-31055 Toulouse, France
Institut de Mathématiques Toulouse, UMR5219 CNRS, 31062 Toulouse, France
François Bachoc
francois.
Instifut de Mathémaliques de Toulouse, UMR5219 CNRS, 31062 Toulouse, France
Jérôme Morio
jerome.morio@onene
ONERA/DTIS, Université de Toulouse, F-31055 Toulouse, France
Timothé Krauth
timothe. @zhaw..r
ONERA/DTIS, Université de Toulouse, F-31055 Toulouse, France
Zurich University of Applied Sciences, Centre for Aviation, Winterthur, Switzerland
Reviewed on OpenReview: https: // openrevieww net/ fo forum? id= n#G9KGssSe
Abstract
Adaptive importance sampling is a well-known family of algorithms for density approxima- tion, generation and Monte Carlo integration including rare event estimation. The main common denominator of this family of algorithms is to perform density estimation with weighted samples at cach iteration. However, the classical existing methods to do so, such as kernel smoothing or approximation by a Gaussian distribution, suffer from the curse of dimensionality and/or a lack of flexibility. Both are limitations in high dimension and when we do not have any prior knowledge on the form of the target distribution, such as its number of modes. Variational autoencoders are probabilistic tools able to represent with fidelity high-dimensional data in lower dimensional space. They constitute parametric family of distributions robust faced to the dimension and since they based on deep neu- ral networks, they are flexible enough to be considered as norpparaitrric models. In this paper, we propose to use variational autoencoder as the auxiliary importance sampling distribution by extending the existing framework to weighted samples. We integrate the proposed procedure existing adaptive importance sampling algorithms and we illustrate its practical interest on diverse examples.
1 Introduction
Importance sampling is a well-known uncertainty quantification method which requires to deal with weighted samples, l.c. a set of observations paired with a set of corresponding weights. estimating an expectation, such as failure probability, in the aim of reducing the variance of the Monte Carlo estimator (Kahn & Harris 1951 Rubinstein & Kroese 2004 Kurtz & Song 2013), for generating points from target probability distribution (Cappé et al. 2004 Cornuet et al.) 2012 Marin et al. 2019). What these algorithms have in common is that they all require to estimate a target probability distribution with weighted samples, and obviously, the accuracy of the algorithm depends on the quality of the imation of the distribution (Chatterjee & Diaconis] 2018). first way for this estimation to use non- non-paran non-parametric models, such as kernel smoothing (Wand Jones 1994 Scott & Sain (005). These models are flexible but they strongly suffer from the curse of dimensionality. Another solution is to use parametric families of distributions, such the Gaussian (Rubinstein & Kroese] 2004) Gaussian mixture (Kurtz & Song.
classically used for

---

Published in Transactions on Machine Learning Research (01/2024)
Geyer et al., 2019 ones, which are more robust in medium-high dimension. However, they sometimes require some prior knowledge on the target distribution, and their lack of flexibility and the huge number of parameters to estimate (Au & Beck 2003) can negatively impact the quality of the estimation when the dimension is high.
2013
In order combine both flexibility and robustness to the dimension, we suggest. here use as the auxiliary importance sampling distribution distribution parameterised by a variational autoencoder, whose main principle has been introduced in the last decade (Kingma & Welling 2014 Kingma et al., 2019). Varia- tional autoencoders are deep generative models for approximating high-dimensional complex distributions of observed data and generating new samples. The specific. feature of a variational autoencoder compared to other density estimation methods is that it performs a dimensionality reduction into a lower dimensional latent space in order to facilitate the estimation. Moreover, in opposition to other dimensionality reduction techniques such as principal component analysis (Wold et al, 1987) or autoencoders (McClelland et al., 1987 Bank et al., 2023), variational autoencoders have good generation properties and give explicitly the approximating distribution, allowing to perform Monte Carlo simulations. This tool is now popular in the machine leaning community but not much in uncertainty quantification.
In the present article, we extend the existing framework to the case of weighted samples by introducing a new objective function. The flexibility of the obtained family of distributions makes it as expressive as a non-parametric model, and despite the very high number of parameters to estimate, this family is much more efficient in high dimension than the classical Gaussian or Gaussian mixture families. Moreover, in order to add flexibility to the model and be able learn multimodal distributions, we consider learnable prior distribution for the variational autoencoder latent variables. We also introduce new pre-training procedure for the variational autoencoder to find good starting weights of the neural networks to prevent as much possible the posterior collapse phenomenon to happen. At last, we explicit how the resulting distribution can be combined with importance sampling, and we exploit the proposed procedure in existing adaptive importance sampling algorithms to draw points from a target distribution and to estimate a rare event probability in high dimension on two multimodal problems.
The remainder of this paper is organized as follows. First, Section 2 provides background on importance sampling. Section formally presents the problem of density estimation and provides review the principle of variational autoencoders and some improvements. Then, Section 4 introducs our suggested extension of variational autoencoders to weighted samples. In addition, Section 5 presents the posterior collapse phenomenon and a procedure to handle it. Section 5 iillustrates the practical interest of the proposed procedure on the generation from some target distributions and on the estimation of failure probabilities on two multimodal problems. Finally, Section concludes the present article and gives future research perspectives stemming from it.
2 Background on importance sampling
In this section, we recall the main principle of importance sampling as well as some important requirements. Then, we explicit the main challenges of importance sampling high dimension and some existing remedies to make the method more accurate. At last, we make brief review of existing methods using neural networks for Monte Carlo simulations and adaptive importance sampling.
2.1
General presentation of importance sampling
Importance sampling (IS) is a classical variance-reduction technique which was introduced in Kahn &
Harris
1951
Hammersley & Handscomb, 1961) and massively used for rare event est estimation (Shinozuaa
1986 Papaioannou et al. 2019) in particular. In the case of the estimation of an expectation 1 where v X' C 12 I. is a black-box function, consists rewriting the expectation according to an auxiliary density g X - R4 as Eg (vb (X) w (X)), where w (x) = f(x)/9(x) is the likelihood ratio. To get unbiased estimate, the support of g must contain the support of E X +
1983
Harbitz
2

---

Published in Transactions on Machine Learning Research (01/2024)
s f(x.. The corresponding IS estimator is then given by:
N N
v ((()) w (x
(1)
a.N
n=1
where (x(*) nE [1.N] is an i..d. sample distributed according to the IS auxiliary distribution g. It is consistent and unbiased, and it has zero-variance if and only if g Sopt with Vx E X, 9opt (x) 10 fx) (Bucklew 2004), the condition that &p is non-negative. However, this optimal density cannot be used in practice because the normalizing constant is 7, which is the quantity to estimate. Nevertheless, one can still expect a significative variance reduction by approximating the optimal density by a sub-optimal auxiliary distribution close to Sopt in a certain sense (more dotails will be given in Section 2.2 Morcover, importance sampling can also be used in a Bayesian context for generating samples according to a target distribution known up to a constant (Cappé et al., 2004 Cornuet et al., 2012 Marin et al. 2019 Martino et al.) 2015; El-Laaham & Bugallo 2019; 2021). We refer to (Bugallo et al. 2017 Elvira & Martino 2021 for a more exhaustive presentation of classical importance sampling algorithms.
The common denominator of every importance sampling procedure is that we need to be able to not only sample from the built auxiliary distribution, but also to have access to its PDF values. In the specific case of the estimation of an expectation, it is crucial to get them to be able to compute the importance sampling estimator in equation [| Moreover, in adaptive importance sampling algorithms, the likelihood ratios are roquired at the end of an iteration to be able to perform the next onc.
2.2 High- digeeniimesion importance sampling
As briefly said in the previous sub-section, the main challenge of every importance sampling method is to approximate target distribution by an auxiliary one. As first way to do so is to use non non-parametric
approximating distributions. Various methods based on kernel smoothing with (Zhang) 1996 Morio 2011
2012) and without (Ang el al. 1992 West 1993 Wand & Jones] 1994) weighted samples have been investi- gated. They are flexible enough to be adapted to a large range of problems, especially multimodal problems. However, they strongly suffer from the curse of dimensionality since the size of the required sample to have a good approximation of the target distribution exponentially grows with the dimension and despite some improvements (Silverman 1986: Zhang et al., 2006: Perrin et al. 2018), these models are very inaccurate in high dimension.
Second, a target distribution can also be approximated within a parametric family of distributions. The most classical way to choose the best candidate is to minimise the Kullback-Leibler divergence (Kullback & Leibler 1951 between the target distribution and the candidate one w.r.t. the paramcters of the family. This choice relevant since a relation between the value of the Kullback-Leiiber divergence and the variance of the importance sampling estimator has been established in (Chatterjee & Diaconis 2018).
The most common and convenient parametric auxiliary families are the Gaussian family (Rubinstein & Kroese 2004 De Boer et al. 2005) and the Gaussian mixlure family (Kurtz & Song 2013 Geyer et al. 2019) for multimodal target distribution. However, the amount of parameters to estimate for these distributions makes them inaccurate when the dimension increases, as well as when there is a lack of a prior knowledge on the form of the target distribution such as its number of modes for example. An inaccurate approximation of the target distribution can lead to the weight degeneracy phenomenon (Rubinstein Glynn| 2009), when a few number of likelihood ratios dominates the other ones, and which strongly badly affects the performances of any importance sampling algorithm. Numerical techniques have been investigated in order to make the density approximation more robust faced to the dimension, such as methods to control the likelihood ratios [onidess 2008 Koblents & Míguez 2015; Martino et al. 2018). shrinkage of the covariance matrix (Ledoit & Wolf 2004 Lahham et al. 20019) or hods to reduce the ellective dimension of the problem (Rubinstein & Glynn 2009: Uribe et al. 2021 El Masri et al., (202]). A final method is to consider families of distributions more adapted to high dimension. In the specific case where the input distribution is the standard Gaussian (vMFNM) family of distributions Mardia & Jupp/ 2000 Wang & Song 2016: Papaioannou et al. 2019 is well-suited and more robust. when the dimension is high. However,
distribution, the von Fishee--iii------a
3

---

Published in Transactions on Machine Learning Research (01/2024)
when considering mixtures of vMFNM distributions for multimodal target distributions, the corresponding algorithms require the knowledge of the number of modes.
2.3 Monte Carlo simulations and importance sampling with neural networks
The use of neural networks for Monte Carlo analyses has also already been investigated in the literature. First, the authors of @Papadrakakii et al. 1996 Hurtado & Alvarez 2001 Papadrakakis & Lagaros 2002)
apply neural networks for Monte Carlo simulation and risk analysis. Morcover, the authors of (Nacssoth et al.
Le et al.) 2018 Lawson et al., 2018 use Aatent-variable models as approximating models for sequential Monte Carlo by maximising a lower bound of the log-likelihood to be more efficient in high dimension. Then, the use of normalising flows for high-dimensional Monte Carlo integration and sampling (Müller et al.) 2019 Gao et al. 2020 Arbel et al.] 2021) have been proposed. In addition, the idea of importance sampling has also been considered to propose more accurate objective function for some latent- variall models (Burda
2018
ct al.
2016).
More related to our work, the authors of (Dowling et al.
combine adaptive importance sampling and variational inference whereas the authors of (Wang et al. 2019 propose an adaptive importance sampling
2018
algorithm supported by variational autocncoder constant. At last, the authors of (Dasgupla Johnson 2024) introduce a new algorithm coupling normalizing flows and importance sampling for the estimation of rare event probabilities in high dimension and without any prior knowledge on the form of the failure domain.
sample from a target distribution known up to a
3 Probability density estimation supported by variational autoencoders
Probability density function (PDF) estimation is a major topic of interest in statistics, as the PDF fully characterises the distribution of a continuous random vector. Given an observed dataset (X()) neji.y " of points from the inpul space X C Rd of dimension 2 1, it consists in the construction of an estimate of the true underlying unknown PDF g* : XX + R., There exist several strategies do so. First, parametric methods consist in picking the best reprosentative of the target density 9 within parametric family of densities. A few classical families (Kotz et al. 2004) are often used and can accurately model fair number of phenomena, but their lack of flexibility and/or the high amount of parameters to estimate when the dimension d increases is a limitation. Second, an important class of non-parametric methods is provided by kernel smoothing (Wand & Jones] 1994 Scott & Sain 2005). Despite the larger ficxibility of the corre- sponding approximating models and some improvements Silverman 1986 Zhang et al. 2006; Perrin et al. 2018), kernel smoothing strongly suffers from the curse of dimensionality. One can also mention functional decomposition methods (Wasserman) 2006) or non- non-parametric Bayesian methods (Ghosal & van der Vaart, 2017), but they also suffer from the curse of dimensionality.
A recent strategy to perform probability density cstimation in high dimension is offered by variational autoencoders (VAEs). The principle of this probabilistic tool was first introduced in (Kingma & Welling 2014). We refer to (Kingma et al. 2019) for an exhaustive and up-to-date presentation. Let us then present in this section the general principle of VAEs.
3.1 Latent variable and Bayesian variational inference
Bayesian variational inference FFx & Roberts 2012 consists in approximating the target distribution g" by paramotric model go. where 0 denotes its parameters. In order to facilitate the density estimation, Bayesian variational inference introduces an unobserved latent variable which lies in a lower-dimensional space Z C Rd, with d, < d. Together, both variables define a joint distribution on X x Z. By marginalising over the latent variable z, the distribution on X' is given by:
go (x) = /
go (x, z) dz.
(2)
Nevertheless, the marginalisation over an unknown and unobserved latent variable makes the computation of the integral intractable, and consequence also the direct computation of the estimated distribution
4

---

Published in Transactions on Machine Learning Research (01/2024)
go. However, thanks to Bayes theorem, it is possible to write:
go (z [x)x p(z) 90 (x (z)
90 (x)
where p is the prior distribution on z and where go (x (z) is the likelihood. Being able to infer the true posterior go (z (x) allows then to compute go (x) from equation 3 To do so, we approximate the true posterior by a variational posterior distribution qá (z xx) chosen within an expressive parametric family of distributions paramcterized by Φ
The choice of the dimension of the latent space dz has huge influence on the accuracy of the density approximation and should be the best trade-off between dimensionality reduction and loss of information. Indeed, performing the variational inference in a lower dimonsional subspace Z instcad of in the high- dimensional input space X reduces the number of parameters of qo (z lx) to estimate and makes the process more accurate. However, the dimension of the latent space must. be large enough to correctly catch the structure of the data, and thus the structure of the true underlying distribution g".
3.2 General principle of VAEs
In order to perform the density estimation, VAE consists in using two neural networks, probabilistic encoder denoted Es parameterized by the weights s and a probabilistic decoder denoted Da parameterized by the weights 0, to model respectively the variational posterior distribution and the likelihood function. Hence, in the case of VAEs, the parameters & and 0 discussed in Section 3.1 are the weights of Es and Do.
First, the probabilistic encoder performs the dimensionality reduction described above. Indeed, given an input point from the imput space x E X, it returns the parameters of the approximating variational parametric distribution ds (z x). Classically, a Gaussian distribution with diagonal covariance matrix is considered as the variational posterior distribution, such that q6 (. (x) N (ux, Ex); where its parameters (μxx Σx) E.g. (x) are the output of the encoder.
Second, the probabilistic decoder sets the parameters of the likelihood function go (x (2). Indeed, it maps point from the latent space E Z into the corresponding parameters of the likelihood. Once more, for continuous data, the most classical choice for the likelihood is the Gaussian distribution with diagonal covariance matrix. Thus, the likelihood distribution is given by go (. [z) ~ V (uz Ez), where its parameters (Hz, E,) = Do (z) are the output of the decoder. As a result, according to equation 2 and with go (x, z) p (z) go (x [z), the estimated distribution go can be seen as an infinite mixture of Gaussian distributions.
Even if the distribution estimated by VAE belongs theoretically to a parametric family, its flexibility, its ability to approximate complex distributions and its form make it closer to non-parametric model. As discussed at the beginning of Section 3 other existing paramartit models strongly face the curse of dimensionality whercas, as shown in the numerical results in Scction D VAEs are much more efficient in high dimension.
The main challenge now is train both neural networks. The most common method for this is to minimize the Kullback-Leibler divergence (Kullback & Leibler 1951) DKL (g"ll between the target distribution g" and the candidate one ge w.r.t. the parameters 0. One can show that it is equivalent to maximise the
expectation Ey. [log (go (X)) This is also called cross-cntropy minimisation (Rubinstein & Kroesc 2004) Nevertheless, the computation of this expectation requires the computation of the integral in equation 2 which is intractable. However, for fixed x E X, one can find a more convenient and easier to compute lower bound of the log- likeeiilioood using the latent variable z (Kingma & Welling 2014 Kingma et al. 2019
log (go (x)) 2 (x)) [log (go (x Z)) DKL (4p (Z|x) |lp (Z))
(4)
Taking the expectation of x w.r.t. g" leads to the objective function of the VAE given by:
Ey [log (ge (X)) 2 Eg [E [log (ge [| Ey- [DKL (gs (Z |p (Z)) = ELBO (Φ. 0)
(5)
The training procedure of the VAE aims then to maximize the ELBO (Evidence lower bound) function
according to the weights (Φ, 0). In practice, all the terms of this objective function are estimated using the
5

---

Published in Transactions on Machine Learning Research (01/2024)
nE[LIN] distributed according to 9". More precisely, the ELBO function combines two opposite phenomena, and optimising it consists in finding the best trade-off between them.
observed sample (x(()
First, maximising the reconstruction term ( [log (go (x|Z))] in order to have an accurate re- construction. More precisely, the encoded distributions do (.) associated to each point from the dataset must be separated enough from each others in the latent space such that the decoder is able to reconstruct the correct point X!n) from it.
Second, minimising the regularisation term, which is the Kullback-Leibler divergence DKL (go (Z|x) [pp (Z)) between the prior and the variational posterior distribution over the dataset. Indeed, the decoder is trained from samples distributed d (.Ix()) for all n E [1,N]] .J. Then, when new X samples are generated, the decoder is applied to new z samples distributed according to the prior p. If this prior is too far away from all the ta (+[X(")), the new generated points may not be representative of the true underlying distribution g", because the decoder is applied outside of its training regime.
3.3 Choice of the prior and VampPrior
Let us discuss the choice of the prior distribution on z. Classically, the most common choice for the prior is the standard Normal distribution in dimension dz. This simple choice practically convenient, especially because it gives an analytical expression of the Kullback-Leibler divergence DKL. (yo (Z |Pp (Z)) since ds (z)x) is also Gaussian. However, it is not optimal in general, for multimodal problems for example, and can lead to over-regularisation and consequently to poor density estimation performances. A more flexible prior is most of the time required.
Theoretically, by rewriting the ELBO function with two regularisation terms and by maximising it w.r.t. the prior p (Makhzani et al. 2015 Hoffman & Johnson 2016), the optimal prior analytically given by:
D" (=(2)= / 2o (zlx) g" (x) Eg lap (z|X)].
This optimal prior distribution called the aggregated posterion. One can note that the aggregated pos- terior is the continuous mixture of all the variational posterior distributions over the whole input space. In particular, it depends on the parameters of the encoder &. A natural way approximate p" is pamp (z) 1 E 9 (z)X(").. Nevertheless, since N is usually large, this empirical distribution is in- tractable in practice and can also lead to over-fitting. To overcome these issues, there are several types of priors to approximate the optimal one p² (Dilokthanakul et al., 2016 Nalisnick et al. 2016 Kim & Mnih 2018 Chen et al. 2018 Lavda et al., 2019 Takahashi et al., 2019 Kalatzis et al. 2020 Negri et al. 2022). each them with advantages and drawbacks. Here, we use the Variational Mizture of Posteriors prior, or VampPrior, introduced in (Tomezak Welling 2018 It consists approximating the optimal prior by a mixture distribution of the form:
K
Pu Pui... ux & (z) = K Σ a (z)uz),
(7)
where K 2 1 is the number of components of the mixture chosen by the user and where the points (ux)ke(1.K] X are learnable pseudo-inputs from the input space. More precisely, the VampPrior distribu- tions with K components constitute a parametric family of distributions parameterised by (u,...., ux, Φ).
Practically, single neural network denoted VP> and parameterised by the weights À maps the vectors of the canonical basis (et) )e(L.)] of IR* into the pscudo-inputs (ux) Therefore, the VampPrior distribution only parameterised by the weights À and d. This new neural network has to be trained in order find the best pscudo-inputs that maximise the performances of the VAE, and these new paramoters have then to appear in the ELBO function which thus becomes:
EIBBO (g, 0. A) Ky [Exa((Ix) [log (go (x|Z))] Eg [DKL. (gs (2(X) I (Z))),

---

Neverthe- Published in Transactions on Machine Learning Research (01/2024)
1 where PX (z) K R== 96 (zIVP> (e.)). The VampPrior is flexible enough to be adapted to many kinds
of problems and has the advantage to depend on the weights of the encoder Ø, as in the form of the optimal prior in equation O
4 Variational autoencoder with weighted samples
In this section, we present a new procedure to estimate, with VAE, a target probability density g", no longer by using observations drawn from g" itself, but rather by using observations drawn from another probability density f X - R+. To the best of our knowledge, this problem has only been investigated (Wang al., 2019) (see Remark 5.1 and Section for comparisons). To that end, we show how to adapt the VAE framework presented in Section 3 to address this issue.
In the context of the rest of the paper, we assume that we have at our disposal dataset (x AE[I,M] E ENN distributed according to a distribution : x + R+ known up to constant, and that we aim estimate the target probability density g" also known up constant. More precisely, the distributions and g" are given by:
{ 1/00 (x)
/cx g (x) /e,
(9)
where f and g" are fully known non-negative functions and where ef and Co are positive constants. At last, we also assume that the support of f contains the support of g"
In the same way as in Section Bl in order to approximate accurately possible the target distribution g", aim to minimise the Kullback-Leibler divergence DRL (g ([9e) w.r.t. the parameters 0, which is equivalent to maximise Eg- [log (go (X))]. However, we have here sample distributed according to f and it is not possible to compute this expectation under g". Fortunately, thanks to the importance sampling trick Kahn & Harris 1951 , it is possible to rewrite it as an expectation under f:
Eg- [log (go (X)) E g (X) (X)
log (go (x))
(10)
The latter expectation still requires the computation of the intractable integral of ge in equation 2 less, by multiplying both sides of equation 4 by the positive weight g (x) /f (x), we can obtain a computable lower bound of the weighted log- liglliicl depending on the latent variable z:
g" (x) f (x)
log (go (x)) 2 9" (x)
(Boo(.1x) [log (go (x|2))] DKL. (go (Z(x) ||x (Z)))
(11)
f (x)
Finally, applying the expectation w.r.t. X f, we get:
9" (X)
Eg [log (go (X)) = E,
log (go (x))
(x)
g" (X)
(X) (Ey (log (go (X|Z)) - DKI. (go (Z(X) (px (2))
F
(12)
As in Section 3 the training procedure of the VAE will here aim maximise the lower bound of equation 12 over the data. At last, even though f and g" are known only up to positive constants Cy and Cg. one ignore these constants in the expression of the function to maximise. Therefore, the objective function in the current framework is given by the new weighted ELBO function:
WELBO (0.0,A) EE E,
g"
Bxx(, [log (go (X/Z))
E,
g [(x)
DKL (gs (Z(X) IPx (21)
(13)
Maximising wELBO (Φ, A) w.r.t. the parameters Φ, 0 and is equivalent to maximising the right-hand side of equation 12 w.r.t. Φ and X, which is lower bound of Eg [log (go (X)))- To sum up, only
7

---

Published in Transactions on Machine Learning Research (01/2024)
by modifying the objective function of the classical VAE from Section g it is possible to learn the high- dimensional target distribution 9" from samples distributed according to an initial distribution f with a VAF. In practice, all the terms of this objective function are estimated using the observed sample (x(")) (n)) nEIL.NI
distributed according to f.
5 New initialisation procedure of the weights of the neural networks to prevent posterior collapse
5.1 Posterior collapse
A classical problem that badly affects the accuracy of the density estimation and the generating properties of the VAE is posterior collapse. As described in many articles (Bowman et al. 2015 Higgins et al. 2016 Sonderby et al.] 2016; He et al. 2019), posterior collapse generally refers to an regulaaguaaiiiation of the VAF, a loss of information, which mathermatically traduces that the Kullback-Leibler term of the objective function vanishes, that DKL (9ø (Zix) ||px. (Z)) 0 for every
X. In order words, every variational posterior distribution collapses into the prior. Getting stuck in a local maxima during the optimisation can
explain this phenomenon (Sonderby et al. 2016 Lucas et al. 2019). Notably, trained VAE affected by posterior collapse is not able to catch the different modes of a multimodal distribution. Note that (Takida ett al. 2020 investigates an alternative definition of posterior collapse based on the mutual information.
The most classical solution to deal with posterior collapse is to introduce a coefficient 3 E [0, 1] before the Kullback-Leibler term of the objective function in order to reduce the regularisation effect (Bowman et al. Higgins al., 2016; Sønderby et al., 2016), but the corresponding objective function no longer a lower bound of the log-likelihood of the model. Others papers propose alternative formulation of the objective function (Rezende & Viola, 2018 Alemi et al.) 2018) whereas the works of (Razavi et al. 2019 Xu & Durrett 2018) investigate new choices for the family of distributions for the prior or the variational
2015
posterior.
5.2 New pre-initialisation procedure
In the present article, we introduce new solution to handle posterior collapse. As explained above, posterior collapse is caused by getting stuck in a local maximum of the objective function. In order to prevent it, we propose a pre-initialisation procedure of the weights of the neural networks Φ, 0 and À in order to start the training of the VAF from more adapted starting points (b("), g(0) and (
First, we initialise the weights by supervised learning. To do so, we randomly pick without replace- ment a sequence of indices (s(k))x=juxI of integers of [[1, N] with probabilities proportional to the family
(5 (x(") /f (x())) nE[L,N] in order to create the sub-sample (x(+(()) kE[1,R] Then, we train the neural
network VPA such that it maps each vector into the corresponding picked point from the dataset (s(k)) by minimising the mean square error:
(VPx (ex) - x(O()) T (VPx (ex) ((()))
(0) = arg min
(14)
The proposed sampling strategy of the sub-sample of size K ensures that it is distributed according to the target distribution g" rather than the initial distribution f. Therefore, once pre-trained, the initial pseudo-inputs (0) VP (ex) will already be representative of the target distribution g", As result, the characteristics of the target distribution g", such as its modes, will be alrcady caught by VP
Second, we initialise the weights & and 0 by unsupervised learning. To do so, in order to already ensure good reconstruction properties, we train the pair encoder soooder/dene /ecoder (Eg, De) as a classical autoencoder on the whole dataset by minimising the mcan square crror between the data points and the reconstructed ones (encoded and then decoded). Once more, thanks to the importance sampling trick, we are able rewrite the loss function as expectation under Moreover, since classical autoencoder only cares about the
8

---

Published in Transactions on Machine Learning Research (01/2024)
mean values Mx and Mz, add penalisation term in order to also pre-train the weights of the networks corresponding to the scale parameters of the encoder in order to initialise the diagonal terms of Σx close to 1 (for normalized data).
We write (E% (x), (x)) = Es (x) and (DB (z) DE (z)) = Do (z) and we recall that E (x) and DØ (z) are diagonal matrices of size respectively d₂ x dz and d d representing the output covariance matrices. We also let log" (D) be obtained by taking the square logarithm of a diagonal matrix D, and we get:
(0 O(O))
arg min Eg. (x-
D/ E (x))) (x- (x-D% DD
E (x)(x)
1 Ir (log" (E (x))
dg
Φ.0
= arg min E 6.0
9" (X)
((x- D' EB (x))) (x- D2 ( (x)(x)))
+
1
Ir (log² (E (x))
5 (X)
g" (X) (X)
((x- Di B (EX=D((E(x))) (x- D/ + (x)
1 a. Tr (log" (o (x))
arg min E,
+
(15)
0.0
Minimising the penalisation term is equivalent then to make each component of the diagonal of E (X) close to 1 for every point of the data set. Since we pre-train the pair encoder /deodrr as a classical autoencoder, we decided to add this penalisation term for normalised data to make sure that this part of the network is also pre-trained.
At last, the full proposed training procedure of the VAE with weighted samples is given in Algorithm Note that this procedure can easily be adapted to the classical VAE framework described in Section 3 only by removing the likelihood ratio g" (x) /f
Algorithm VAE-IS
Require: (x()) nE[L,N] J,K,dz
1: Randomly pick K points (x()(A)) kE[l,K] within the dataset according to the procedure described at the
beginming of Section 5.2
Train VP> by minimising equation
14 and get A (0)
3: Train (Eş, De) by minimising equation
15 and get (d) ((0)
Train the whole VAE (E,, Do. VPx) by maximising equation starting from (6 (())) (0)
return Trained VAE.
Remark 5.1. The objective of [Wag et al. 2019) is the same as ours: learning a target distribution g" by using observations from an initial distribution f with a VAE. The main principle of their algorithm the following:
1. Irain a classical VAE as in Section 3 with the available observations from f,
2. draw new points from the trained VAE and compute the corresponding likelihood ratios,
3. boostrap among the new sample with probabilities proportional to the likelihood ratios.
They iteratively repeat this procedure and the final sumple is theoretically distributed according to distri- bution close to g" Contrary to their work, we propose one-step procedure to achieve the same goal by modifying the objective function of the VAE. Moreover, it does not seem clear how to have access to the PDF of the resulting approzimating distribution of g" in (Wang et al. [201]] (if at all possible) because they only have access lo a boolstrapped sample, whereas il is struighlforward with our approach. In addition, we also introduce a new pre-training procedure of the weights of the neural networks in order to prevent posterior collapse. At last, they only test their algorithm on low dimensional (10 and 6) uni-modal target distributions whereas, as developed in Section 0 we test our procedure on more complex test cases. In Section 6 our procedure is compared with the algoritlom of Wang et al. (201)), and is more efficient.
9

---

Published in Transactions on Machine Learning Research (01/2024)
6 Application to high-dimensional non-parametric importance sampling
In this section, we first show how to apply the proposed procedure to the classical framework of importance sampling. Then, in order to illustrate the practical interest of the previous efforts, we evaluate numerically the performances of the suggested procedure and we compare them to the performances of some existing methods. The code reproduce the numerical experiments publicly available at: Jullom6661//mportaace- Saaptpprrraaace-nmiinnngg
com/
Sampling- VAE.gis
6.1 Importance sampling supported by a VAE
In this paper, we propose to select the IS auxiliary distribution within the distributions parameterised by VAEs. In the IS framework, the target distribution corresponding 9 of Section here 9opt, which is known up to constant. Morcover, we assume that the initial sampling distribution is perfectly known. Then, in order compute the likelihood ratios w of the estimator in equation let explain how to get the corresponding PDF of the resulting distribution go- The most naive way to do so is (Wang el al. 2019):
1. draw sample (x(m), (()) "E[L.N] according to the joint distribution go (x,z) z) p() go (x (z)
X x Z,
2. estimate the PDF values associated to the sample with:
N
g0
(((m)
g0 (x(m) (n) z(1)
(16)
However, the use of these estimated values of the PDF makes the IS estimator equation 1 biased since go (X()) is the denominator. Morcover, the PDF values (go (x(")) nE[1,]] are not independent. Thus, the estimator in equation D is no longer a sum of independent random variables and as a result, its convergence is no longer guaranteed.
Hence, in order keep the convenient properties of the classical IS estimator in equation (unbiasodness and convergence), we adopt the following procedure:
approximate the whole marginal distribution go in equation by:
90 M (.)
1
Σ go Z (m)
(17)
M
me[L.M] is an i.i.d. sample from the latent space Z distributod according to the prior distribution p.
where M 2 and where m))
2. draw sample (X (x())) NEIINN according to gM and compute the estimator og NN
The main point of this new procedure is that M for which we can generate new points and compute exactly and straightforwardly the PDF for all x E X. 90
longer work with go at all but only with the distribution
In that case, we come back to the classical and well-known IS framework from Section and g.M allows to keep the convenient statistical properties of the IS estimator.
Proposition 6.1. For every fized M 2 1, the IS estimator FIS is unbiased and convergent as N +00.
96 N
Proof. See Appendix
The detailed architecture of the VAE used for the following test cases shown in Figure 1 each test case, we K 75 and M 103, and we will explicit the dimension d, of the latent space chosen for each example.
Morcover, for
10

---

Published in Transactions on Machine Learning Research (01/2024)
Pseudo-inputs
(e kE[LA]
Hx
VPa
Pa.
kx
Fx
H4
(x )EELNNI f
Dg
90 ~ g
Σx
Σz
Latent Space
Figure Reprosentation of our suggested VAE architecture for IS in high dimension.
6.2 Generation of samples by adaptive IS
A first application of density estimation with weighted samples is adaptive IS for sampling from target distribution g² Sometimes, sampling directly from the target distribution g" is not possible, especially when it is only known up to constant (Ghosal & van der Vaart 2017). Then, starting from imitial proposal distribution f. adaptive IS algorithms consist in iteratively approximating the target distribution within a family of proposal distributions and at the end, the resulting sample is expected to be drawn according to a distribution close to g"
In both examples, we use classical adaptive IS scheme (AIS) as (Marin al. 2019), except that we return the weighted sample generated only at the last iteration. We use distribution parameterised by VAE as the proposal distribution each iteration (AES.VAE). The detailed pseudo-code of the procedure is described in Algorithm 2
Algorithm AIS-VAL VAE
Require: f. N, K, dz, M, ilmax
1: Draw
x ())
according to f
nE[1.N]
2: Compute wn g (x) /s x
for il E [1, mmax] do
Train VAE (E Do, VPx) according to the procedure in Algorithm [ with the weighted sample (x(it-)) wn 1))
Compute the estimate of the resulting distribution 90. M
6: Draw X(1r
(it)
according to
[I,1] 7: Compute wn = g" (x /99...
( (it)
end for
return Weighted sample (x((i)
(it) ((t)
11

---

Published in Transactions on Machine Learning Research (01/2024)
6.2.1 Example in dimension 10
First, let us consider a 10-dimensional bimodal target distribution given by:
si x (2. x I10 +N (-2.5 x 11,,110)
(18)
where 110 = (I,..., 1)² E R and Io is the 10-dimensional identity matrix. For comparison purposes, we execute as well the AIS algorithm with a Gaussian mixture with two components (AIS-GM) as the proposal distribution. We also execute the algorithm introduced in Wang et al. (2019) with the proposed setting (AIS- Wang ct al. (2019) ) as well as the algorithm introduced in this article but with standard Gaussian prior and without the pre-training procedure of Section 5.2 (ATS-standard VAE). We will also use the Kullback- Leibler divergence (Kullback & Leibler 1951] to quantitatively compare the performances of every method, since it is a classical measure to quantify the dissimilarity between two distributions. However, its value strongly depends on both compared distributions and also on their dimension. Therefore, it is not possible to have a clear idea of the expected scale and it is only relevant to compare multiple values. The standard Gaussian distribution in dimension 10 is the starting proposal distribution for both algorithms. We perform Prep 100 executions of each algorithm with 10 iterations. Then, at each iteration, we draw N 10 new points, and the dimension of the latent space chosen for this problem is =4
For every algorithm, the results are as follows: either both modes, and so the whole target distribution, are perfectly found and approximated, or only one mode is found but it is well approximated. Examples of the representation of the final approximating sample are given in Figure 2 However, the main difference between the four algorithms is their success rate, i.c. the frequency with which the resulting distribution finds both modes. In the first place, as shown in Table [ the AIS-VAE algorithm finds both modes more than 70% of the time and gives very good approximation of the target distribution. In the second place, the AIS-GM algorithm only finds both modes less than 10% of the time over the 'erp 100 repetitions. Third, the algorithm of Wang et al. (2019) is intermediate because it finds both modes 46% of the time. These results show that a distribution parameterised by a VAE is more likely to catch the characteristics of a target distribution than a single Gaussian or a Gaussian mixture, and even though the latter knew the number of modes in advance. Morcover, the use of a tempered version of both the encoder and the decoder as in Wang et al. (2019) definitely helps to explore larger area of both the latent space and the input space and thus to catch the characteristies of target distribution, compared to the AIS procodure with standard VAE with a standard Gaussian prior and without the pre-training procedure as the IS auxiliary sampling distribution. The latter has a success rate of 25%. Nevertheless, the use of a standard Gaussian prior prevents to reach better performances on a multimodal problem as our full procedure. Thus, these numerical results strongly highlight the interest of the use of both VampPrior and of the pre-training procedure for multimodal problems.
The likelihood ratios correspouding to the samples adaptively guide the algorithm from the initial distribution f to the target one gî. However, even more so when the dimension increases, one or some values can be significantly larger than the others, and as result the adaptive IS algorithm will go in that direction. Nevertheless, we can not explain why this phenomenon seemss such less important when using a VAE rather than a Gaussian mixture. A first possibility is that the dimensionality reduction performed by the VAE makes the process more robust. A second one can be the differences between the two learning procedures themselves. The learning procedure is based on mini-batches for the VAE whereas the whole dataset is used in single step for learning the Gaussian mixture.
Table 1: Comparison of the four algorithms. The first row of the table represents the success rate of the corresponding algorithm and the quantity (DKL) (DKL)success ysuccess represents the mean value of the Kullback-Leibler divergence between the empirical generated distribution and the target of over the successful samples.
AIS-VAE
AIS-GM
AIS Wang et al. 2019)
ATS-standard VAE
Success rate
72%
7%
46%
25%
(DKL.) mean
youreess
2.48 x 10-2
1.13 X 10
4.69 x 10
1.21 X 10
12

---

Published in Transactions on Machine Learning Research (01/2024)
33-
02
0 BE25
76
X1
X1e
X7
Ki
(a) Example for the AIS- VAE algorithm with 2 modes found.
(b) Example for the AIS- VAE algorithm with
mode found.
X
X1
X
X
34
02
(c) Example for the AIS-GM algorithm with modes found.
(d) Example for the A1S-GM algorithm with 1 mode found.
Figure 2: Graphical representation as 10 histograms of the 10 marginals of the final generated sample for each case for the VIS-VAE and ATS-GM algorithms. The continuous orange line represents the PDF of the marginals the target distribution gi. The continuous green line represents PDF of the marginals of the starting distribution f.
6.2.2 Example dimension 20
Second, let consider 20-dimensional target distribution 9i defined by the marginal distributions and the dependence structure given in Table 2 The starting distribution chosen here is f f= = (02o, 2 x I20). We perform 10 iterations and we draw N = 10 new points at each one. The dimension of the latent space for this problem is d2 = 8. Examples of the resulting empirical distribution returned by the AIS-VAE algorithm can be found in Figure 3 Moreover, since the target distribution g has only one mode, is interesting note that the performances of the AIS-standard VAE algorithm (not reported in Figure 3 are exactly the same the ones of the full proposed procedure in Algorithm [ We can see graphically that the target distribution is most of the time well approximated by the resulting distribution. However, because of very high negative values of the log-PDF, it is not possible to compute numerically the Kullback- Leibler divergence between the empirical distribution and the target distribution 9i, and then to use it as a quantitative measure of the quality of the approximation. Indeed, we recall that the expression of the forward Kullback-Leiblee divergence between the target distribution 9i and the approximating one Ô2 is DKL (gāllg) = Eg; [log (gž (X)) - log (ō2 (X))]. For estimating this quantity, we first draw an i..d. sample (X(1) ne]I.N] from L2, compute the log-PDF values for each generated point and we return the empirical mean. However, if for some E [] the PDF Ō2 (X())) is very very close to O, the corresponding log-value log (ōz (x ())) will be very large negative value. Then, it might be possible that the software, Python for us, not handle such high value and consider it "Not number' (NaN). If so, the estimated Kullbac--Leibler divergence will also be NaN.
Moreover, for comparison purposes, we applied the AIS algorithm with a single Gaussian (AIS-SG) as the proposal distribution, since of has only one mode. However, does not work at all with the same setting, Indoed, the starting distribution algorithm. Then, it necessary choose of closer to g2, by modifying the mean vector for example. Once more, this example shows that using a distribution parameterised by VAE, allows to have larger flexibility
N (020, 2 I20) is too far away from the target Ø for the AIS-SG
13

---

Published in Transactions on Machine Learning Research (01/2024)
C 2
RRCER SRRCE
(a)
(b)
4
20
(c)
(d)
Figure 3: Graphical representation as 20 histograms of the 20 marginals of four realisations of the final generated sample for the AIS- VAE algorithm. The continuous orange line represents the PDF of the marginals of the target distribution gi. The continuous green line represents the PDF of the marginals the starting distribution f.
14

---

Published in Transactions on Machine Learning Research (01/2024)
Table 2: Student (v, μ4, o): I-dimensional Student distribution with V > 0 degrees of freedom, with mean μ E and scale paramcter a O. LogN a): distribution of exp (4) with A a 1-dimensional Gaussian random variable of mean μ E R and scale a > O. Triangular (a, m, b): 1-dimensional Triangular distribution where a b are the lower and upper bounds and where E [a, b] is the mode. The correlation matrix is given by Bag 1 (i- j) + 1/4 x 1 (li jl = 1) for all (i.j) E [1,20]
Input
Distribution
X1
Student (4, -2.
X2
LogN (0, 1)
X3
Triangular (1,3, 5)
X1 to X20
.V (2,
Copula
Normal copula with correlation matrix R E R 20020
than with a classical single Gaussian distribution, because the corresponding adaptive IS algorithms have a larger range. At last, we also applied the algorithm proposed in Wang et al.] (2019). However, it does not work at all neither. Indeed, as explained in Section 2 at each iteration of an adaptive importance sampling algorithm, we need to compute the likelihood ratios. We recall that the expression to compute the PDF values of the resulting distribution in their algorithm is given in equation [G Then, it seems that this expression a sum makes the computation of the PDF values highly numerically unstable when the dimension increases. Some of them are infinite and as result, an error arises.
Finally, we point out that because of the weight degeneracy phenomenon (Rubinstein & Glynn 2009) which occurs in IS when the dimension increases, the AIS-VAE algorithm, in the same manner as the other adaptive IS algorithms, will not work very well in higher dimension than 20.
6.3 Adaptive importance sampling for estimating failure probabilities
Another application ion of density estimation with weighted samples is adaptive IS for reliability analysis. Re- liability analysis consists in the estimation of the failure probability pt = TP (v (X) > t) = Ey [1 (vb (X) > t)] for a fixed known threshold & E R. Classical Monte Carlo sampling is not adapted to this problem when Pt is getting smaller because its computational cost becomes too large to obtain an accurate estimation. Therefore, several techniques have been developed in order to reduce the variance of the estimation of pe: one can mention FORM/SORM methods (Hasofer & Lind 1974 Breitung 1984), subset sampling (Cérou et al. 2012) or line sampling (Koutsourelakis et al. 2004) for example.
In the IS framework, because of the rarity of the failure event, it is challenging to directly pick the best representative of 9opt (x) x 1 (bb (x) > t) f (x) within the parametric family of densities by minimising the KL divergence to 9ope- To overcome this issue, the authors of (Rubinstein Kroese [004 De Boer et al., 2005) introduced the so-called multi-level cross entropy method. At each iteration of the algorithm, the cross-entropy problem is solved for less rare intermediate failure event until reaching the true onc. The intermediate events are set from quantile parameter p. Several parametric families of densities have been used for this algorithm in the literature: the Gaussian family (Rubinstein & Kroese 2004 De Boer et al. 2005), the Gaussian mixture family (Kurtz & Song 2013 Geyer et al. 2019) and the vMFNM family (Wang & &Song 2016 Papaioannou et al. 2019) for example. For the following two test cases, we use the classical multi-level cross-entropy but with auxiliary densities parameterised by VAFs (CF-VAF) and in both of them, we set the dimension of the latent space to d. = The detailed pscudo-code of the procodure is described in Algorithum 3
At each iteration of the algorithm, we draw No. 10a points and the quantile parameter is set to p= 0.25 for the first example and to p = 0.15 for the second. We perform nrep = 100 realisations of cach estimator to have an estimation of the error. First, we compare the performances of Algorithu X with the cross-entropy algorithm using a VAE with a standard Gaussian prior and without the pre-training procedure of Section as the IS auxiliary distribution (CE-standard VAE). Sccond, since the Gaussian and the Gaussian mixture families do not perform well at all in high dimension, we will compare the CE- VAE algorithm with the multi-level cross-entropy algorithm with mixtures of vMFNM as auxiliary densities (CE-vMFNM). As said
15

---

Published in Transactions on Machine Learning Research (01/2024)
Algorithm 3 Cross-entropy algorithm with VAE as the auxiliary density
Require: N, p 00,1, K, dz, M, it =0
1: Draw
x nE[I,N]
according to
2: Evaluate (1) = v (0)
x (0) (n) for all n [1,
Sort 1(1) - ... v V(n) and set min (1 (((I-)N):.))
4: Compute wn (0) 1
yn (0)
70)
5: while a <t do
6: it = itt
Train VAE (Eø, Do, VP>) according to the procodure in Algorithm |
with the weighted sample
(it-1) Wh (()-11))
Compute the estimate of the resulting distribution go M according to 90.t
9: Draw
X
(i))
NE[LNN]
10: Evaluate yn = vb
for all n E [1,N]
11: Sort (it)
(at) (1) S ... S ((N) and set Vit = min (1((--)
12:
Compute w = 1 Un
(it)
Yet f (x (it)
00.0
M
[ (ir)
13:
end while
14: return N-1 DN E w
in Section 2.2 the latter algorithm requires the number components of the mixture, so we will test different setups. For comparison purposes, let define the following quantities:
No the number calls to the function for each execution of the algorithm,
mmean; the mean estimated failure probability over the nrep rcalisations,
COV (fr): the coefficient of variation of the estimator over the Trep realisations,
PMC: a coefficient. allowing to compare the efficiency of the method compared to the classical Monte Carlo method defined by PMC NMC /Not where NMC - p.) the required sample to reach the same coefficient of variation COV (p.) with the classical Monte Carlo method. If PMC 1, the method is more efficient than Monte Carlo and conversely.
(p, (p:)* (1--:)/pp..o
the size of
6.3.1 Four branches
First, let. us consider the analytical problem
CChiron el al.
2023
given for any even dimension d2 2 by:
1
Di=
34
d
Σd D4
Vx Rd, v1 (x) min
Vd
(19)
>4 d/2 23; - >i
2//d/2+1 Di
Vd
-1' 24 +
Ei 24
Vd
The input vector associated to this function the standard Gaussian distribution in dimension The failure threshold is set 3.5 such that the theoretical value the failure probability p. 9.3 10-4, This problem is challenging because has failure modes. graphical representation of this function dimension d = 2 well the failure limit state represented in Figure 4a Here, we set the dimension
16

---

Published in Transactions on Machine Learning Research (01/2024)
2
x
2
3,5
-4
-10
-10
-5
10
-4
2:
X1
(b) Representation of the two-dimensional PDF of the prior distribution in the latent space associated to the last iteration of one execution of the CE- VAE algorithm.
(a) Two-dimensional representation of the function t as well as the failure limit state corresponding to 3.5 in red.
Figure
of the problem to 100, and the numerical performances of the algorithms are given in Table 3 By comparing both the coefficient of variation and the coefficient PMC associated to each estimator, we see that the proposed CE- VAE algorithm provides better performances than the CE-vMFNM algorithm here, and without any prior knowledge on the number of failure modes. Figure 4b represents the PDF of the prior distribution in the latent space at the last iteration for one execution of the CE-VAE algorithm. One can clearly see that the four failure modes are well represented it. This a clear illustration of the interest of both VampPrior and of the pre-training procedure on a high-dimensional multimodal problem. Therefore when generating new points the four modes will be represented. Moreover, one can note that the choice of the number of components in the mixture for the CE-+MFNM algorithm is crucial. Indeed, since the current problem has failure modes, need at least 4 components in the mixture in order to get good results. The algorithm for 3 components does not always converge because it sometimes reaches the maximal number of iterations, and the estimation error is large. In addition, unreported visual inspection shows that the empirical distribution of the Nrep = 100 realisations of the CE-vMFNMA algorithm has a heavy tail, and therefore explains the observed error. At last, the results of the CE-standard VAE algorithm are a clear illustration of the huge negative impact of the posterior collapse phenomenon on the performances of the algorithm. Indeed, without VampPrior and the pre-training procedure, the standard VAE never finds any of the four modes and as result, the corresponding algorithm never converges since it always reaches the maximal number of iterations.
Table 3: Comparison of the CE- VAE algorithm with the CE-vMNFM algorithm with different numbers of components and the CE-standard VAE algorithm for the four branches problem. The number after the "CExMMFNM" acronym represents the number components of the mixture given as an algorithm input.
CE-VAE
CE-MMENMB
CE-MFENA
CE-MMFNM5
CE-standard VAE
Ntot
40000
88000
50000
50000
200000
Smean
9.310 X 10
1.319 x 10
9.835 x 10
9.315 x 10
9.446 x 10
COV (p,)
5.31%
512.8%
31.3%
7.56%
34.83%
PMC
9.54
4.64 x 10-4
2.19 x 10
3.76
4.43 x 10-2 2
17

---

Published in Transactions on Machine Learning Research (01/2024)
6.3.2 Duffing oscillator
The second example is the Duffing oscillator introduced in (Zuev]
2009];
and we consider here its writing
discretised in the Trequency domain as in (Shinozuka & isodat
1991
Papaioannou et al. 2019). The
function of interest is the maximal displacement of the oscillator at amax
2sec given by:
2 (x) min {ur (Amax), u (Amax) u2]
(20)
where #n 0.1 and U2 -0.06, and where the displacement u (a) of the oscillator satisfies for all a 2 0: mü (a) +
d/2
+ cu (a) + &o (u(o) + ta((aa) -mo
Σ
(x, cos (wra) + Xay2++ sin (w,a)),
(21)
1=1
with m 1000 kg, c == 200 Ns/m, = 1000 (2T)²² N/m, y= m- wi iAw, Aw = 30m /d and
V0.01AW. Moreover, the initial conditions are set to (0) and à (0) 1.5, the random variable X follows standard Gaussian distribution of dimension d 200 and the failure threshold is set to f = 0. It
is possible to verify that there are then two failure modes. The reference value of the failure probability is
Pref = 4.28 10-4 and it has been computed with large Monte Carlo estimation of size 10 (Chiron et al. 2023
Table 4: Comparison of the CE-VAB algorithm with the CE-vMNFM algorithm with different number of components and the CE-standard VAE algorithm on the Duffing oscillator problem. The number after the "CE-MMFNM" acronym represents the number of components of the mixture given as an algorithm input.
CE-VAE
CF-vMENMI
CF-vMFNM2
CE-MMENM3 CR-standard VAE
Ntot
30000
200000
40000
40000
47900
Pt
omicn
4.27 x 10
4.09 x 10
4.25 x 10
4.26 x 10
3.39 x 10 4
COV (pc)
8.69%
45.85%
3.72%
2.51%
31.09%
PMC
10.3
5.55 x 10
42.2
92.8
5.04 x 10
The numerical performances of the algorithms are given in Table 4 By comparing once more both the coefficient of variation and the coefficient VNC associated to each estimator, the CE-VAE algorithm provides satisfying performances, but the CE-vMFNM2 and CE-vMFNM3 algorithms outperform it. Indeed, their coefficient of variation between two and three times smaller than the one of the CE-VAE algorithm and although the CE-vMFNM algorithms require one more iteration converge, their associated coefficient vNC is much better. The current problem is probably very well suited for this family of auxiliary distributions and can explain these very good results in such high dimension. However, already said in the previous example, contrary to the proposed CE-VAL VAE algorithm, the CE-vMFNM algorithm requires the knowledge of the number of components in the mixture, which can often be a limitation. Here, at least two components are required because we can observe that the CE-vMFNMI algorithm does not work at all since it reaches the maximal number of iteration for each execution. Then, the proposed training procedure of the VAE allows identify both failure modes without any prior knowledge and to get satisfying results anyway. At last, in the same way the previous example, the CE-standard VAE algorithm often suffers from the posterior collapse phenomenon and does not perform well as the CE-VAE algorithm.
7 Conclusion
In the present article, we interested in probability density estimation with weighted samples, i.e. the estimation of a target distribution 9" by using a sample drawn from another distribution f. It is major topic of interest in isaiis which has many applications, some of them presented in Section 0 We suggest to approach g" by a distribution parameterised by a VAE. To do so, we extend the existing VAE framework to the case of weighted samples, by introducing the new objective function wELBO maximise in order to learn the best parameters of the neural networks. Even if the corresponding distribution theoretically belongs to parametric. family, its characteristies make it closer to non-parametric model. Despite the
18

---

Published in Transactions on Machine Learning Research (01/2024)
very high number of parameters to estimate, this family is much more efficient in high dimension than the classical Gaussian or Gaussian mixture families. Moreover, in order to add flexibility to the model, and more precisely to be able to learn multimodal distributions, we use a learnable prior distribution for the latent variable called VampPrior. We also introduce a new pre-training procodure for the VAE in order to find good starting points (o ( ) for the maximisation of wELBO and to prevent as much as possible the posterior collapse phenomenon to happen. At last, we illustrate and discuss the practical interest of the proposed method. We first describe the classical IS framework for estimating an expectation and we show how to exploit the resulting distribution to that purpose. Then, we apply the proposed procedure in an adaptive TS algorithm for drawing points according to a target distribution. Both examples in dimension 10 and 20 show that a VAE has a larger flexibility and is more likely to catch the specificity of a distribution than a Gaussian or a Gaussian mixture distribution. Finally, we introduce the proposed method in an adaptive IS algorithm for estimating a failure probability in high dimension. We observe that the resulting estimation is always quite accurate without any prior knowledge on the form of the problem, in opposition to the vMNFM family. These numerical examples also strongly illustrate and support the use of VampPrior and of the pre-training procedure for multimodal problems.
It is important to note that the computational time of the training of neural networks, and so of a VAE, is growing with the number of training points and the input dimension. As comparison, the training time of the other parametric families used in this article, not relying on neural networks, is almost instantaneous whereas one execution of the proposed procedure can last from less than one minute to minutes on a CPU. However, once the VAE is trained, the generation of new points is instantancous, and when considering time-expensive black-box functions, the training time of a VAE becomes negligible. Moreover, the suggested suggest here a unique architecture quite simple
architecture in Figure I is not the only possibility. Indeed, to implement but general enough to be adapted for every considered problem. However, if it is necessary, an interested user can modify the suggested architecture to be more adapted to its data.
Finally, this article shows that a VAE is a powerful tool for density estimation with or without weighted samples, and we apply it for adaptive IS. Then, one can investigate the use of VAEs for others algorithuss or uncertainty quantification problems requiring high-dimensional density estimation, such as an improved adaptive IS algorithm for failure probability estimation (Papaioannou et al. 2019) or some MCMC methods
CChi & Greenberg
1995
Au & Beck 2001) for example.
References
Alexander Alemi, Ben Poole, Ian Fischer, Joshua Dillon, Rif A Saurous, and Kevin Murphy. Fixing broken ELBO. In International conference on machine learning, pp. 159-168. PMLR, 2018.
George L Ang, Alfredo H-S Aug, and Wilson H Tang. Optimal importance-sampling density estimator. Journal of enginoering mechanics, 118(6):1146-1163, 1992.
Michael Arbel, Alex Matthews, and Arnaud Doucet. Annealed flow transport Monte Carlo. In International Conference on Machine Learning, pp. 318 330. PMLR, 2021.
Siu-Kui Au and James L Beck. Estimation of small failure probabilities in high dimensions by subset simulation. Probabilistic engincering mechanics, 16(4):263-27) 2001.
Siu-Kui Au and JL Bock. Important sampling in high dimensions. Structural safety, (5 2003.
Dor Bank, Noam Koenigstcin, and Raja Giryes. Autocncoders. Machine Learning for Data Science Handbook: Data Mining and Knowledge Discovery Handbook, pp. 353 374, 2023.
Samuel R Bowman, Luke Vilnis, Oriol Vinyals, Andrew M Dai, Rafal Jozefowicz, and Samy Bengio. Gener- ating sentences from a contimuous space. arXiv preprint arXiv:1511.06349, 2015.
Karl Breitung. Asymptotic approximations for multinormal integrals. Journal of Engineering Mechanics, 110(3):357-300 366, 1984.
James Bucklew. Introduction to rare event simulation. Springer Science & Business Media, 2004.
19

---

Published in Transactions on Machine Learning Research (01/2024)
Monica F Bugallo, Victor Elvira, Luca Martino, David Luengo, Joaquin Miguez, and Petar M Djuric. Adaptive importance sampling: The past, the present, and the future. IEFE Signal Processing Magazine, 34(4):60 79, 2017.
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov. Importance weighted autoencoders. Proceding in International Conference on Learning Representation, 2016.
Olivier Cappé, Arnaud Guillin, Jean-Michel Marin, and Christian P Robert. Population Monte Carlo, Journal of Computational and Graphical Statistics, 13(4):907-929, 2004.
Frédéric Cérou, Pierre Del Moral, Teddy Furon, and Arnaud Guyader. Sequential Monte Carlo for rare event estimation. Statistics and Computing, 22(3):795 908, 2012.
Sourav Chatterjee and Persi Diaconis. The sample size required in importance sampling. The Annals of Applied Probability, 28(2):1099 1135, 2018.
Ricky TQ Chen, Xuechen Li, Roger B Grosse, and David K Duvenaud. Isolating sources of disentanglement in variational autoencoders. Advances in neural information processing systems, 31, 2018.
Siddhartha Chib and Edward Greenberg. Understanding the Pottropoiis-assting algorithm. The american statistician, 49(4):327-38) 335, 1995.
Marie Chiron, Christian Genest, Jérôme Morio, and Sylvain Dubreuil. Failure probability estimation through high-dimensional elliptical distribution modeling with multiple importance sampling. Reliability Engineering & System Safety, 235:109238, 2023.
Jean-Marie Cornuet, Jean-Michel Marin, Antonietta Mira, and Christian P Robert. Adaptive multiple
importance sampling. Scandinavian Journal of Statistics, 39(4))798-8)) 812, 2012.
Agnimitra Dasgupta and Erik A Johnson. REIN: Reliability Estimation via Importance sampling with Normalizing flows. Reliability Engincering & System Safety, 242:109729, 2024.
Pieter-Tjerk De Boer, Dirk P Kroese, Shie Mannor, and Reuven Y Rubinstein. A tutorial on the cross-entropy method. Annals of operations research, 134(1):19 67, 2005.
Nat Dilokthanakul, Pedro AM Mediano, Marta Garnelo, Matthew CH Lee, Hugh Salimbeni, Kai Arulku- maran, and Murray Shanahan. Deep unsupervised clustering with Gaussian mixture variational autoen- coders. arXiv preprint arXiv:1611.026A 2016.
Matthew Dowling, Josue Nassar, Petar M Djurić, and Mónica Bugallo. Improved adaptive importance sampling based on variational inference. Tn 2018 26th European Signal Processing Conference (EUSIPCO), pp. 1632 1636. IEEE, 2018
Yousef El-Laham and Mónica Bugallo. Stochastic gradient population Monte Carlo. Processing Letters, 27:46 50, 2019.
IEFF Signal
Yousef El-Laham and Mónica F' Bugallo. Policy gradient importance sampling for bayesian inference. JEEE Transactions on Signal Processing, 114424-42566 2021.
Yousef Fl-Laham, Victor Elvira, and Mónica Bugallo. Recursive shrinkage covariance learning in adap- tive importance sampling. In 2019 IEEE 8th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pp. 624 628. IEEE, 2019.
Maxime El Masri, Jérôme Morio, and Florian Simatos. Improvement of the cross-entropy method in high dimension for failure probability estimation through a one-dimensional projection without gradient esti- mation. Reliability Engineering & System Safety, 216:107991, 2021.
Victor Elvira and Luca Martino. Advances in importance sampling. arXiv preprint arXiv:2102.05407 2021.
Charles W Fox and Stephen J Roberts, tutorial on variational Bayesian inference. Artificial intelligence review, 38:85 95, 2012.
20

---

Journal of the Published in Transactions on Machine Learning Research (01/2024)
Christina Gao, Joshua Isaacson, and Claudius Krause. i-flow: High-dimensional Integration and Sampling with Normalizing Flows. Machine Learning: Science and Technology. (44):045023, 2020.
Sebastian Geyer, Jason Papaioannou, and Daniel Straub. Cross entropy-based importance sampling using
Gaussian densities revisited. Structural Safety, 76:15 27, 2019.
Subhashis Ghosal and Aad van der Vaart. Fundamentals of nonparametric Bayesian inference, volume 44. Cambridge University Press, 2017.
John Michael Hammersley and David Christopher Handscomb. Monte Carlo Methods. In Proceedings. volume 7, pp. 17. US Army Research Office., 1961.
Alf Harbitz. An efficient sampling method for probability of failure calculation. Structural safety, 3(2): 109 115, 1986.
Abrahan M Hasofer and Niels C Lind. Exact and invariant second-moment code format. Engineering Mechanics division, 100(1):111-21,, 1974.
Junxian He, Daniel Spokoyny, Graham Neubig, and Taylor Berg-Kirkpatrick. Lagging inference networks and posterior collapse in variational autoencoders. arXiv preprint arXiv:1901. 10055449 2019.
Trina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mo- hamed, and Alexander Lerchner. beta-VAE; Learning basic visual concepts with constrained variational framework. In International conference on learning representations, 2016.
Matthew D Hoffman and Matthew J Johnson. ELBO surgery: yet another way to carve up the variational evidence lower bound. In Workshop in Advances in Approximate Bayesian Inference, NIPS, volume 1, 2016.
Jorge E Hurtado and Diego A Alvarez. Neural-network-bascl reliability analysis: comparative study. Computer methods in applied mechanics and engineering, (1-(1)):113 132, 2001.
Edward L Tonides. Truncated importance sampling. Journal of Computational and Graphical Statistics, 17 02):295-3114 2008.
Herman Kahn and Theodore E Harris. Estimation of particle transmission by random sampling. National Bureau of Standards applied mathematics series, 12:27-30, 1951.
Dimitris Kalatzis, David Eklund, Georgios Arvanitidis, and Søren Hauberg. Variational autocncoders with Riemannian Brownian motion priors. arXiv preprint arXiv:2002.05227, 2020.
Hyunjik Kim and Audriy Mnih. Disentangling by lactorising. International Conference on Machine Learning, pp. 2649 2658. PMLR, 2018.
Diederik P. Kingma and Max Welling. Auto-Eua Iato-Encoding Variational Bayes. Article accepted in the 2nd International Conference Learning Representations 2014, 2014.
Diederik P Kingma, Max Welling, et al. An introduction to variational autoencodors. Foundations and Trends@ in Machine Learning, 12(4):307 392, 2019.
Eugenia Koblents and Joaquín Míguez. A population Monte Carlo scheme with transformed weights and its application to stochastic kinetic models. Statistics and Computing, 25(2):407-425, 2015.
Samuel Kotz, Narayanaswany Balakrishnan, and
Norman L Continuous multivariate distributions, Volume 1: Models and applications, volume 1. John Wiley
Johnson.
& Sons, 2004.
Phadeon-Stelice Stelios Koutsourclakis, HJ Pradlwarter, and GI Schuäller. Reliability of structures in high dimen- sions, part I: algorithms and applications. Probabilistic Engineering Mechanics, 19(4):409-411 417, 2004.
21

---

Published in Transactions on Machine Learning Research (01/2024)
Solomon Kullback and Richard A Leibler. On information and sufficiency. The annals of mathematical statistics, 22(1):79 86, 1951.
Nolan Kurtz and Junho Song. Cross-entropy- adaptive importance sampling using Gaussian mixture. Structural Safety, 42:35 44, 2013.
Frantzeska Lavda, Magda Gregorová, and Alexandros Kalousis. Improving VAE generations of multimodal data through data-dependent conditional priors. arXiv preprint arXiv:1911.1088, 2019.
Dieterich Lawson, George Tucker, Christian A Naesseth, Chris Maddison, Ryan P Adams, and Yee Whye Teh. Twisted variational soquential Monte Carlo. In Third workshop on Bayesian Doep Learning (NeurIPS), 2018.
Tuan Anh Le, Maximilian Igl, Tom Rainforth, Tom Jin, and Frank Wood. Auto-encoding sequential Monte Carlo. Proceding in International Conference Learning Representation, 2018.
Olivier Ledoit and Michael Wolf. A well-conditioned estimator for large-dimensional covariance matrices. Journal of multivariate analysis, 88(2):365-411 2004.
James Lucas, George Tucker, Roger Grosse, and Mohammad Norouzi. Understanding posterior collapse in generative latent variable models. 2019.
Alireza Makhzani, Jonathon Shlens, Navdeep Jaitly, lan Goodfellow, and Brendan Frey. Adversarial autoen- coders. arXiv preprint arXiv:1511.05644, 2015.
Kanti V Mardia and Peter Jupp. Directional statisties, volume 2. Wiley Online Library, 2000.
Jean-Michel Marin, Pierre Pudlo, and Mohammed Sedki. Consistency of adaptive importance sampling and recycling schemes. Bernoulli Society for Mathematical Statistics and Probability, 25(3):1977 1998, 2019.
Luca Martino, Victor Elvira, David Luengo, and Jukka Corander. An adaptive population importance sampler: Learning from uncertainty. IEEE Transactions on Signal Processing, 63(16):4422 4437, 2015.
Luca Martino, Victor Elvira, J Míguez, A Artés- Rodréggzd,, and PM Djurié. A comparison of clipping strategies for importance sampling. In 2018 IEEE Statistical Signal Processing Workshop (SSP), pp. 562. IEEE, 2018.
James L McClelland, David E Rumelhart, PDP Research Group, et al. Parallel distributed processing, volume Explorations the microstructture of cognition: Psychological and biological models. MIT press, 2, 1987.
Jérôme Morio.
Non-parametric adaptive importance sampling for the probability estimation of a launcher impact position. Reliability engineering & system safety, 96(1):178-183, 2011.
Jérôme Morio. Extreme quantile estimation with nonparametric adaptive importance sampling. Simulation Modelling Practice and Theory, 27:76-89, 2012.
Thomas Müller, Brian McWilliams, Fabrice Rousselle, Markus Gross, and Jan Novák. Neural importance sampling, ACM Transactions on Graphics (ToG), 38(5):1 19, 2019.
Christian Naesseth, Scott Linderman, Rajesh Ranganath, and David Blei. Variational sequential Monte Carlo. In International conference on artificial intelligence and statistics, pp. 968-977. PMLR, 2018.
Eric Nalisnick, Lars Hertel, and Padhraic Smyth. Approximate inference for deep latent Gaussian mixtures. In NIPS Workshop on Bayesian Deep Learning, volume 2, pp. 131, 2016.
Marcello Massimo Negri, Vincent Fortuin, and Jan Stuchmer. Meta-leare learning richer priors for VAEs. In Fourth Symposium Advances in Approximate Bayesian Inference, 2022.
22

---

Published in Transactions on Machine Learning Research (01/2024)
Manolis Papadrakakis and Nikos D Lagaros. Deliability-baaeee structural optimization using neural networks and Monte Carlo simulation. Computer methods in applied mechanics and engineering, (32):3344911333 -5007 2002.
Manolis Papadrakakis, Vissarion Papadopoulos, and Nikos D Lagaros. Structural reliability analyis of elastic- plastic structures using neural networks and Monte Carlo simulation. Computer methods in applicd mechanics and engineering, 136(1-2):145-1633) 1996.
Tason Papaioannou, Sebastian Geyer, and Daniel Straub. Improved cross roy importance sampling with a flexible mixture model. Reliability Engineering & System Safety, 191:106564, 2019.
Guillaume Perrin, Christian Soize, and Noura Ouhbi. Data-driven kernel representations for sampling with an unknown block dependence structure under correlation constraints. Computational Statistics & Data Analysis, 119:139-155 154, 2018.
Ali Razavi, Aäron den Oord, Ben Poole, and Oriol Vinyals. Preventing posterior collapse with delta- VAEs. arXiv preprint arXiv:1901 004166 2019.
Danilo Jimenez Rezende and Fabio Viola. Taming VAEs. arviv preprint arXiv:1810.0059, 2018.
Reuven Rubinstein and Peter W Glynn. How deal with the curse of dimensionality of likelihood ratios in Monte Carlo simulation. Stochastic Models, 25(4):547 568, 2009.
Reuven Y Rubinstein and Dirk P Kroese. The cross-entropy method: a unified approach to combinatorial optimization, Monte Carlo simulation, and machine learning. Springer, 133, 2004.
David W Scott and Stephan R. Sain. Multidimensional density estimation. Handbook statistics, 24: 229-261, 2005.
Masanobu Shinozuka. Basic Analysis of Structural Safety. Journal of Structural Enginccring-asoe, 109: 721 740, 1983.
Masanobu Shinozuka and George Deodatis. Simulation of stochastic processes by spectral representation. 1991.
Bernard W Silverman, Density estimation for statistics and data analysis, volume 26. CRC press, 1986.
Casper Kaae Sonderby, Tapani Raiko, Lars Maaløe, Seren Kaae Sonderby, and Ole Winther. Ladder varia- tional autoencoders. Advances in neural information processing systems, 29, 2016.
Hiroshi Takahashi, Tomoharu Iwata, Yuki Yamanaka, Masanori Yamada, and Satoshi Yagi. Variational autoencoder with implicit optimal priors. Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pp. 1066-5073 2019.
Yuhta Takida, Wei-Hsiang Liao, Toshimitsu Uesaka, Shusuke Takahashi, and Yuki Mitsufuji. AR-ELBO: Preventing Posterior Collapse Induced by Oversmoothing in Gaussian VAE. 2020.
Jakub Tomezak and Max Welling. VAE with a VampPrior. In International Conference on Artificial Intelligence and Statistios, pp. 1214 1223. PMLR, 2018.
Felipe Uribe, Jason Papaioannou, Youssef M Marzouk, and Daniel Straub. Cross-entropy- importance sampling with failure-informed dimension reduction for rare event simulation. STAM/ASA Journal Uncertainty Quantification, 9(2):818 847, 2021.
Matt P Wand and M Chris Jones. Kernel smoothing, CRC press, 1994.
Hechuan Wang, Mónica F Bugallo, and Petar M Djuríć. Adaptive importance sampling supported by variational auto-encoder. In 2019 IEEE 8th International Workshop on Computational Advances Multi-Sensor Adaptive Processing (CAMSAP), pp- 619 623. TEFE, 2019.
23

---

Published in Transactions on Machine Learning Research (01/2024)
Ziqi Wang and Junho Song. Cross-entropy-basee adaptive importance sampling using Mises-Fishar Fisher mixture for high dimensional reliability analysis. Structural Safety, 59:42-52, 52, 2016.
Larry Wasserman. All of nomparametric statistics. Springer Science & Business Media, 2006.
Mike West. Approximating posterior distributions by mixtures. Journal of the Royal Statistical Socicty:
Series (Methodological), 55(2):409 422, 1993.
Svante Wold, Kim Esbensen, and Paul Geladi. Principal component analysis. Chemometrics and intelligent laboratory systems, 2(1-3):37 52, 1987.
Jiacheng Xu and Greg Durrett. Spherical latent spaces for stable variational autoencoders. ar preprint arXiv:1808.1088 2018.
Ping Zhang. Nonparamctric importance sampling. Journal of the American Statistical Association, 91(45)): 1245 1253, 1996.
Xibin Zhang, Maxwell King, and Rob J Hyndman. A Bayesian approach to bandwidth selection for multivariate kernel density estimation. Computational Statistics & Data Analysis, 50(11):3009 3031, 2006.
Konstantin Zuev. Advanced stochastic simulation methods for solving high-dimensional reliability problems. Hong Kong University of Science and Technology (Hong Kong), 2009.
A Proof of Proposition 6.1
Let M 2 1 and let FM be the o-algebra generated by the i.i.d. random variables Z(.,..., Z(M) distributed according to the prior p. We recall that g/ (.) M- Σm= 9o (.Iz(m)) and that the expression of 7 is given in equation T
9.N
A.1
Unbiasedness
We have:
E
p 90
= Ep
E
fIS 90
FM
.N
f (x(1))
Eyp
N Σ (1)
Fas
D6dd (X(1))
N
f (x(())
Ep
E
(x())
Fas
by linearity of the conditional expectation
D0 (x(1n)
M
(X)
Ep
E
M
9ø M
(X)
FM
conditionally to FM,
(n)
an i.idd. sample distributed according to gø
ne[1.N]
= Ep Ey (bb (x)) by importance sampling
Ep [7]
Then, 718 26 N is an unbiased estimator of I.
24

---

Published in Transactions on Machine Learning Research (01/2024)
A.2 Convergence
By the law of large numbers, conditionally to FMr, we have:
f (x(())
7
(n))
a.s.
I.
(22)
9B
.N
N
9o (Xin)) N --0
n=1
Then,
P
s 96
1
js " iN
.N
N
(1
Ep
M
sgr NN No0000 I
Fas
Epp
P
7 % ,N N
FM
Ep [1]
1.
So the estimator 7
%g .N converges almost surely to I.
25


Reviews:
[{"review_id": "EA0AsTBYn4", "reviewer": "Reviewer_G7sn", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\nOverall\n- The claims of the paper matches the contents very well.\n- The contents of the paper are not only useful for people interested in the particular setting it studies (learning density using weighted samples, IS, failure estimation) but could be a good reading for general audience in probabilistic machine learning to learn how various techniques can be used together in practice.\n\nWriting\n- Each part of the paper is very well-written with sufficient background knowledge introduced. I believe this would be very helpful for people who are not familiar with one of the related filed.\n\nMethod\n- The proposed method is neat and based on a set of new or established techniques that are technically sound.\n\n**Questions**\n\nMethod\n- Unbiasedness: Does (17) actually lead to unbiased IS? It looks to me that it's still biased because of the use of Monte Carlo + g appearing in the denominator. Its improvement from (16) is mostly the independence. Can the authors clarify on this?\n\n**Weaknesses**\n\nWriting\n- Abstract: The abstract reads more like a list of \"things we do\" (or concrete contributions) that should appear in the end of the introduction. For abstract, I would like to see a more high-level summary/storyline of the challenges and intuitions of the proposed, e.g. starting with why and then what. If the authors agree with my view, perhaps they could rewrite that part of the abstract.\n- Overall storyline: The current storyline looks a bit weird to me. In particular, the introduction of VampPrior seems to be not well-motivated, until I see it is needed in later sections. Perhaps the authors should consider bring the technique challenges (e.g. VampPrior is motivated by (17) for indepedence) in section 5 upfront and then motivate individual choices.\n- Related work: Although they in the form they are proposed can not work with weighted samples, normalizing flows should at least be mentioned (e.g. in the beginning of section 2) as they are considered as the best performing method for general density estimation.\nExperiments\n- While I agree that methods based on mixture of Gaussians requires prior knowledge of the number of modes, it's also a valid argument to say one can just pick a larger, conservative number. Maybe that would hurts the performance if the number is too large? If that's true (from new experiments), maybe that's something worth including.\n\n**Misc**\n\nPresentation\n- Consider making table 1 to the top instead of in the middle of the main texts.\n- Figure 4 and 5 are way too large\n- Somehow I have difficulty selecting texts on the PDF file.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper demonstrates how to learn the optimal proposal distribution in importance sampling (IS) using variational autoencoder (VAEs) and use it in adaptive IS. This requires a few key ingredients: \n- First, the learning of VAEs is done using weighted samples, for which the paper introduces a new learning objective called weighted ELBO. \n- Second, the paper uses VampPrior in VAEs such that the construction of density estimators can be effectively dependent on the prior, making the estimates independent as required by IS. \n- Additionally, a new way based on supervised + unsupervised pre-training is introduced to avoid posterior collapse, which is needed as the problems considered as multi-modal for which VAEs can suffer from posterior collapse. \n- Finally, the paper shows how to use these together to perform adaptive IS. \n\nExtensive simulations in multiple setups are performed to illustrate, validate and compare the proposed method against established baselines. Some of the limitations (e.g. still suffering from high-dimensional problems) are also discussed.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Requests\n- *critical* See my comments about the unbiasedness. This needs to be clarified and addressed accordingly.\n- *good-to-have* See my comments about writing in *Weakness* above. Those can strengthen the paper.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "qzlMWo7kdT", "reviewer": "Reviewer_VqJg", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n\nI think the proposed loss function is **very** nifty. Specifically, the prospect of training a generative model using data from a different distribution is very cool and I think could lead to a lot of cool applications.\n\n\n# Weaknesses \n\nWhile I like the proposed approach, the paper has a lot of weaknesses. First, I think the writing can be greatly improved. As the novel loss function relies on the importance sampling trick **and** many of the results presented are a combination of the VAE + AIS, I was surprised that the background section didn't go over IS or AIS. I think this would greatly improve the flow of the paper. Also, there was virtually no related works section. I found this very surprising given that the combination of VAEs and IS-based methods is not new, i.e., [1, 2]. Moreover, there are a slew of methods that try to update the proposal of AIS by minimizing a loss function [3, 4]. The related works section would only strengthen the paper it would help highlight the main contribution which is that IS is done over the data used to train the VAE **and not the latents** which is what has been done in previous works. Also, the formatting instructions for TMLR state to use `\\citep` and `\\citet` but both are missing in the paper. While this is a minor point, it does reduce the readability of the paper.\n\nNext, the weighted ELBO in equation 12 is **not** a lower bound of the expected log marginal likelihood, $E_{g^*}[\\log g_\\theta(x)]$ nor does it lead to the same gradient as the lower bound in eq 11. Although the normalization constants are not a function of the parameters, it does scale the ELBO and need to be taken into account to ensure that you're optimizing a lower-bound to the expected log marginal likelihood; if the normalization constants were instead an additive term then you could throw them away and would still be optimizing the expected log marginal likelihood. Instead, the proposed lower-bound (and the computed gradients) is a lower-bound to a scaled expected log marginal likelihood, $\\frac{c_f^{-1}}{c_g^{-1}} E_{g^*}[\\log g_\\theta(x)]$. Below, I demonstrate this (where $\\mathcal{L}(x)$ is the ELBO, $g^*(x) = \\tilde{g}^*(x) c_g^{-1}$, and, $f(x) = \\tilde{f}(x) c_f^{-1}$)\n$$ \\log g_\\theta(x) \\geq \\mathcal{L}(x) $$\n$$ \\frac{\\tilde{g}^*(x)}{\\tilde{f}(x)} \\log g_\\theta(x) \\geq \\frac{\\tilde{g}^*(x)}{\\tilde{f}(x)} \\mathcal{L}(x) $$\n$$ E_f \\left[ \\frac{\\tilde{g}^*(x)}{\\tilde{f}(x)} \\log g_\\theta(x) \\right] \\geq E_f \\left[ \\frac{\\tilde{g}^*(x)}{\\tilde{f}(x)} \\mathcal{L}(x) \\right] $$\n$$ \\int f(x) \\frac{\\tilde{g}^*(x)}{\\tilde{f}(x)} \\log g_\\theta(x) dx \\geq \\int f(x) \\frac{\\tilde{g}^*(x)}{\\tilde{f}(x)} \\mathcal{L}(x) dx $$\n$$ \\int \\tilde{f}(x) c_f^{-1} \\frac{\\tilde{g}^*(x)}{\\tilde{f}(x)} \\log g_\\theta(x) dx \\geq \\int \\tilde{f}(x) c_f^{-1} \\frac{\\tilde{g}^*(x)}{\\tilde{f}(x)} \\mathcal{L}(x) dx $$\n$$ \\int c_f^{-1} \\tilde{g}^*(x) \\log g_\\theta(x) dx \\geq \\int c_f^{-1} \\tilde{g}^*(x) \\mathcal{L}(x) dx $$\n$$ \\int \\frac{c_f^{-1}}{c_g^{-1}} c_g^{-1} \\tilde{g}^*(x) \\log g_\\theta(x) dx \\geq \\int \\frac{c_f^{-1}}{c_g^{-1}} c_g^{-1} \\tilde{g}^*(x) \\mathcal{L}(x) dx $$\n$$ \\int \\frac{c_f^{-1}}{c_g^{-1}} g^*(x) \\log g_\\theta(x) dx \\geq \\int \\frac{c_f^{-1}}{c_g^{-1}} c_g^{-1} g^*(x) \\mathcal{L}(x) dx $$\n$$ \\frac{c_f^{-1}}{c_g^{-1}} E_{g^*(x)} \\left[ \\log g_\\theta(x) \\right] \\geq \\frac{c_f^{-1}}{c_g^{-1}} E_{g^*(x)} \\left[ \\mathcal{L}(x) \\right] $$\nThis was just glossed over by the proposed lower bound but it is important to highlight, especially since the authors highlight that beta-VAEs don't optimize a lower bound to the expected log marginal likelihood.\n\nIn section 4, the authors introduce a pre-initialization procedure to mitigate posterior collapse, which I think is great to include, but isn't well-motivated. Posterior collapse is a problem, but in the VampPrior paper, the authors showed empirical evidence that using the VampPrior mitigates posterior collapse in hierarchical deep generative models, which are much more prone to posterior collapse compared to standard deep generative models. Thus, did the authors experience posterior collapse and they found this pre-initialization procedure mitigated this? Empirical evidence is needed to demonstrate that this procedure mitigates posterior collapse, or at the very least leads to improvement gains compared to not using the pre-initialization procedure.\n\nIn section 5.1.2 the authors state the procedure used in Wang et al., 2019. leads to a biased IS estimator as the estimated pdf values appear in the denominator; they then propose an alternative that \"retains the convenient statistical properties of the IS estimator\". I don't think this is a fair criticism. While it is true that the approach used by Wang et al., 2019 leads to a biased IS estimator of $\\frac{f(x)}{g_\\theta(x)}$ the proposed approach is also a biased estimator of this ratio. The proposed approach leads to an unbiased estimate of $\\frac{f(x)}{g_\\theta^M(x)}$ but as $g_\\theta^M(x)$ is an estimator of $g_\\theta(x)$, it has the same problem as Wang et al., 2019. While a priori, I expect the proposed approach would lead to lower variance, both are biased.\n\nAlso, it isn't clear how this algorithm is combined with AIS. Is the VAE retrained every AIS iteration using samples from the previous iteration? Is the VAE retrained from scratch or is warm-started? The paper would greatly benefit from an algorithm table.\n\nLastly, my biggest concern is the results section. First, the lack of comparisons is not acceptable. This is especially true given that the proposed approach is very similar to Wang et al., 2019. I think at the minimum, all results should be compared against 1) Wang et al., 2019 2) the proposed method without the initialization procedure, and 3) any AIS algorithm that adapts the proposal by optimizing some form of loss function, i.e. [3]. Second, error bars are missing for the first 2 experiments. Third, in section 5.2.2, the authors state that the KL can't be computed due high values of the pdf; I'm confused by this statement. What is high in values, the target? Is it forward KL or reverse KL? If it is forward KL, then I don't see how the target being high in value (assuming that it is the target) would lead to a numerically unstable computation of KL. Fourth, the combination of proposed VAE with cross-entropy minimization for failure probabilities is non-trivial and I think at the very least needs a section in the appendix. For example, in the duffing oscillator, what is the likelihood function? \n\n# References\n[1] Importance Weighted Autoencoders. https://arxiv.org/abs/1509.00519   \n[2] Variational Sequential Monte Carlo. https://proceedings.mlr.press/v84/naesseth18a/naesseth18a.pdf  \n[3] Stochastic Gradient Population Monte Carlo. https://ieeexplore.ieee.org/ielaam/97/8966529/8903438-aam.pdf   \n[4] Policy Gradient Importance Sampling for Bayesian Inference. https://ieeexplore.ieee.org/abstract/document/9468927", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this work, the authors propose a novel algorithm for building a proposal to do importance sampling (IS), and its adaptive counterpart adaptive importance sampling (AIS). Specifically, the proposal is parameterized by a deep generative model trained using an extension to the VAE framework. Specifically, a novel loss function is proposed that allows a deep generative model to be trained using data from a different distribution.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Below, I list my requested changes\n\n1. Include background info on IS and AIS in the appropriate background section. **Critical**\n2. An explicit related works section. **Critical**\n3. Fix citations. **Minor**\n4. Address the loss function concern I raised above **Critical**\n5. Provide empirical evidence of the usefulness of pre-initialization and posterior collapse **Critical**\n6. Rewrite the biased estimator claim in section 5.1.2 **Critical**\n7. Revamp results section with comparisons **Critical**", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "WSxhMkfsxC", "reviewer": "Reviewer_89qk", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper does a good job of providing a thorough background with clear explanation on VAEs. \n\nIt’s nice that the paper is structured to provide pointers to related work throughout, rather than having a dedicated section to Related Work. \n\nI think this paper could use some more motivation to help explain the problem setup (see requested changes). \n\nI think the authors should spend more time explaining why PDF estimation is an important problem. The authors state in Remark 4.1 that Wang 2019 does not provide a method to compute the PDF of the desired distribution g^*, but it is unclear to the reader why this is important. The authors should remind the reader that they wish to apply their method to the problem of importance sampling. \n\nThe results in section 5.2.1 are interesting, but I find it a weakness of the paper that the phenomenon is not thoroughly understood. I also think there should be some ablations demonstrating which part of the authors approach is necessary to make their method work (e.g. how important is the VampPrior). It would be interesting to see for example if the success of the authors method can be predicted by the result of the pre-training initialisation of the VampPrior and/or the encoder and decoder. An experiment in 2 or 3 dimensions would be helpful, perhaps the phenomenon can be fully explained there.\n\nI also think there should be further explanation for the KL divergence metric used in experiment 5.2.1. This is the only metric that measures the sample quality from the joint distribution. It would be nice to provide some context/interpretation for the metric here, since right now it is hard to interpret its scale (the only takeaway I have right now is that it is lower, hence better, for the authors method compared to the baseline).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper begins by giving a thorough background on VAEs, as well as the Variational Mixture of Posteriors (VampPrior) in Section 2.\n\nSection 3 introduces the methodology and the problem set up. The problem setup is to estimate the normalising constant of an empirical distribution using only samples from another empirical distribution (which is also known up to a constant). The authors propose training a VAE with VampPrior using the importance sampling trick to tackle this problem.\n\nThe authors propose a pre-training setup to avoid the problem of posterior collapse. They pre-train the VampPrior, initialising its outputs to be close to points with large ratio g/f.\n\nThe authors use their method to provide improved importance samplers. They test their improved importance samples on several tasks in statistics.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "problem setup. For example, why is PDF estimation relevant vs. sampling from the distribution? It’s stated briefly in the introduction that the authors apply their approach to important sampling, but I think it would be helpful to continually remind the reader of the problem the authors are trying to solve. \n\nIt would be helpful to motivate why it is important to sample points according to the ratio g/f when pre-training the VampPrior.\n\nIn autoencoder pre-training (eq. 14), can you further explain the loss term for the diagonal covariance matrix? For example, I think 1) It would be helpful to remark that the output of the encoder/decoder is the log variance and 2) that the encoder/decoder only output d_z entries corresponding to the diagonal of the covariance matrix. \n\nAlgorithm 1 should explain how the k points are sampled to pre-train the VampPrior. \n\nPotential typo in the paragraph after equation 15 when stating optimality conditions for g_opt. Should be \\forall x \\in \\mathcal{X}.\n\nIn Figure 1, I think it’s unnecessary (i.e. crosses a level of abstraction we don’t need to) to show the Encoder and Decoder are parameterised as MLPs.  \n\nI found Section 5.1.2 confusing. It would be helpful to expand this section and provide more exposition.\n\nIn Section 5.2.2, can you explain why large values of the PDF prevent computing the KL divergence?", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All three reviewers found that the authors' rebuttal was convincing and it addressed all of their issues. Thus, they unanimously lean towards acceptance of the paper. Overall, the paper is a nice contribution to TMLR and I second the reviewers' sentiment towards acceptance.

Reviewrer #VqJg mentioned some actionable suggestions (mostly for the experimental section) which would be great to include in the final version of the paper.
Decision: Accept (poster)