Particle Semi-Implicit Variational Inference

Authors: jen ning lim
adam m johansen
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-03 13:40:30

Abstract: Semi-implicit variational inference (SIVI) enriches the expressiveness of variational
families by utilizing a kernel and a mixing distribution to hierarchically define the
variational distribution. Existing SIVI methods parameterize the mixing distribution
using implicit distributions, leading to intractable variational densities. As a result,
directly maximizing the evidence lower bound (ELBO) is not possible, so they
resort to one of the following: optimizing bounds on the ELBO, employing costly
inner-loop Markov chain Monte Carlo runs, or solving minimax objectives. In this
paper, we propose a novel method for SIVI called Particle Variational Inference
(PVI) which employs empirical measures to approximate the optimal mixing
distributions characterized as the minimizer of a free energy functional. PVI arises
naturally as a particle approximation of a Euclidean–Wasserstein gradient flow and,
unlike prior works, it directly optimizes the ELBO whilst making no parametric
assumption about the mixing distribution. Our empirical results demonstrate that
PVI performs favourably compared to other SIVI methods across various tasks.
Moreover, we provide a theoretical analysis of the behaviour of the gradient flow
of a related free energy functional: establishing the existence and uniqueness of
solutions as well as propagation of chaos results.

Main paper: 
 introduction: 
 In Bayesian inference, a quantity of vital importance is the posterior p(x|y) = p(x, y)/ p(x, y)dx, where p(x, y) is a probabilistic model; y denotes the observed data, and x the latent variable. An ever-present issue in Bayesian inference is that the posterior is often intractable. This is because the normalizing constant is available only in the form of an integral, and approximation methods are required. One popular method is variational inference (VI) (Jordan, 1999; Wainwright et al., 2008; Blei et al., 2017) . The essence of VI is to approximate the posterior with a member from a variational family Q where each element of Q is a distribution q θ (called "variational distribution") parameterized by θ. These parameters θ are obtained via minimizing a distance or discrepancy (or an approximation of it) between the posterior p(•|y) and the variational distribution q θ .Here, we focus on semi-implicit variational inference (SIVI) (Yin and Zhou, 2018) . It enables a rich variational family by utilizing variational distributions, which we refer to as semi-implicit distributions (SIDs), defined asEQUATIONwhere k : R dx × R dz → R + is a kernel satisfying k(x|z)dx = 1; r ∈ P(R dz ) is the mixing distribution and P(R dz ) denotes the space of distributions with support R dz , with its usual Borel σ-field, with finite second moments: E r [∥z∥ 2 ] < ∞. Here, and throughout, we assume that the distributions and kernels of interest admit densities. SIDs are very flexible (Yin and Zhou, 2018) and can express complex properties, such as skewness, multimodality, and kurtosis. These properties might be present in the posterior but typical variational families may fail to capture them. There are various approaches to parameterizing these variational distributions: current techniques utilize neural networks built on top of existing kernels (e.g., Gaussian kernels) to define more complex kernels (Titsias and Ruiz, 2019) , and/or utilize pushforward distributions (a.k.a., implicit distributions (Huszár, 2017)) (Yin and Zhou, 2018) . On choosing a parameterization, an approximation to the posterior is obtained by minimizing the exclusive Kullback-Leibler (KL) divergence. This optimization has the same solution as minimizing the free energy (or the negative evidence lower bound) defined as E(k, r) := log q k,r (x) p(x, y) q k,r (dx).(2)However, since the integral in q k,r is typically intractable, directly optimizing E is not feasible. As a result, SIVI algorithms focus on designing tractable objectives by using upper bounds of E (Yin and Zhou, 2018) ; expensive Markov Chain Monte Carlo (MCMC) chains to estimate the gradient of E (Titsias and Ruiz, 2019); and optimizing different objectives such as score matching which results in min-max objectives (Yu and Zhang, 2023) .In our work, we propose an alternative parameterization for SIDs: kernels are constructed as before (with parameter space denoted by Θ) whereas the mixing distribution r is obtained by optimizing over the whole space P(R dz ). We motivate the case for minimizing a regularized version of the free energy E denoted by E λ (see Eq. ( 4)); thus, SIVI can be posed as the following optimization problem: arg min (θ,r)∈Θ×P(R dz ) E λ (θ, r). As a means to solving the SIVI problem, we construct a gradient flow that minimizes E λ where the space Θ × P(R dz ) is equipped with the Euclidean-Wasserstein geometry (Jordan et al., 1998; Kuntz et al., 2023) . Via discretization, we obtain a practical algorithm for SIVI called Particle Variational Inference (PVI) that does not rely upon upper bounds of E, MCMC chains, or solving minimax objectives.Our main contributions are as follows:(1) we introduce a Euclidean-Wasserstein gradient flow minimizing E λ as means to perform SIVI;(2) we develop a practical algorithm, PVI, which arises as a discretization of the gradient flow that allows for general mixing distributions;(3) we empirically compare PVI compared with other SIVI approaches across toy and real-world experiments and find that it compares favourably; (4) we study the behaviour of the gradient flow of a related free energy functional to establish existence and uniqueness of solutions (Prop. 8) as well as propagation of chaos results (Prop. 9).The structure of this paper is as follows: in Section 2, we begin with a discussion of previous approaches to parameterizing SIDs and their relationship with one another. Then, in Section 3, we show how PVI is developed: beginning with designing a well-defined loss functional, the construction of the gradient flow, and how to obtain a practical algorithm. In Section 4, we study properties of a related gradient flow; and, in Section 5, we conclude with experiments to demonstrate the efficacy of our proposal. For sake of brevity, we defer our discussion of related works to App. A. 
on implicit mixing distributions in sid: 
 This section outlines existing approaches to parameterizing SIDs with implicit distributions and how these choices affect the resulting variational family. Before we begin, we shall summarize the key assumptions of SIVI. The kernel k is assumed to be a reparametrized distribution in the sense of Salimans and Knowles (2013) ; Kingma and Welling (2014) ; Ruiz et al. (2016) . In other words, the kernel k is defined by the pair (ϕ, p k ) where ϕ : R dz × R dx → R dx and p k ∈ P(R dx ) such that k(•|z) = ϕ(z, •) # p k Furthermore, to ensure that it admits a tractable density, the map ϵ → ϕ(z, ϵ) is assumed to be a diffeomorphism for all z ∈ R dz with its inverse map written as ϕ -1 (z, •). From the change-of-variable formula, its density is given as k(•|z) = p k (ϕ -1 (z, •)) | det ∇ x ϕ -1 (z, •)|. We sometimes write k ϕ,p k to denote the underlying ϕ and p k explicitly. Furthermore, the kernel k is assumed to be computable and differentiable with respect to both arguments.Several approaches to the parameterization of SID have been explored in the literature. One can define the variational family by choosing the kernel and mixing distribution from sets K and R respectively, i.e, the variational family is Q(K, R) := {q k,r : k ∈ K, r ∈ R}. Yin and Zhou (2018) focused on a fixed kernel k with r being a pushforward (or "implicit") distribution, i.e., r ∈ {g # p r : g ∈ G} =: R G;pr where G is a subset of measurable mappings from the sample space of p r to R dz . Thus, the Q YiZ -variational family is Q({k}, R G;pr ). On the other hand, Titsias and Ruiz 
particle variational inference: 
 In this section, we present our proposed method for SIVI, called particle variational inference (PVI). Similar to prior SIVI methods, the algorithm utilizes kernels (denoted by k θ ) with parameters Θ which satisfy the assumptions listed in Section 2. One example is k θ ∈ K Θ;ϕ,p k where Θ is a function space induced by a neural network. We slightly abuse the notation Θ to also indicate its corresponding weight space R d θ . The novelty of this algorithm is that, for the mixing distribution, we directly optimize over the space P(R dz ) which loosens the requirement for the neural network in the kernel to learn complex mappings. The result is a "simpler" optimization procedure and increases expressivity over existing methods. Thus, the variational parameters of PVI are (θ, r) ∈ Θ × P(R dz ) =: M with its corresponding variational distribution defined as q θ,r := k θ (•|z)r(z)dz. PVI arises naturally as a discretization of a gradient flow minimizing a suitably defined free energy on Θ × P(R dz ) endowed with the Euclidean-Wasserstein geometry (Jordan et al., 1998; Ambrosio et al., 2005; Kuntz et al., 2023) . In Section 3.1, we begin by constructing a suitably defined free energy functional; then, in Section 3.2, we formulate its gradient flow; finally, in Section 3.3, we construct PVI from its gradient flow.As with other VI algorithms, we are interested in finding variational parameters that minimize (θ, r) → KL(q θ,r , p(•|y)). This optimization problem can be cast equivalently as: r) , where E : M → R : (θ, r) → q θ,r (x) log q θ,r (x) p(x, y) dx.arg min (θ,r)∈M E(θ,(3)Before we can solve this problem, we must ensure that it is well-posed. In other words, it must admit minimizers in M. In the following proposition, we outline various properties of E: Proposition 2. Assume that the evidence is bounded log p(y) < ∞ and k is bounded; then we have that E is (i) lower bounded, (ii) lower semi-continuous (l.s.c.), and (iii) non-coercive.The proof can be found in App. E.1. Prop. 2 tells us that even though E possesses many of the properties one looks for in a meaningful minimization functional, it lacks coercivity (in the sense of Dal Maso (2012, Definition 1.12)): a sufficient property to establish the existence of solutions. The key to showing non-coercivity is that we can construct a kernel k θ (x|z) that does not depend on z.At first glance, this issue might seem contrived but we note that this problem is closely related to the problem of posterior collapse (Lucas et al., 2019; Wang et al., 2021) . To address non-coercivity, we propose to utilize regularization and define the regularized free energy as:EQUATIONwhere R λ is a regularizer with parameters λ. In Prop. 3, we show that if R λ is sufficiently regular, then the E λ enjoys better properties than its unregularized counterpart E. Proposition 3. Under the assumptions of Prop. 2, if R λ is coercive and l.s.c., then the regularized free energy E λ is (i) lower bounded, (ii) l.s.c., (iii) coercive. Hence it admits at least one minimizer in M.The proof can be found in App. E.2. From here forward, we shall focus on regularizers of the formR E λ : (θ, r) → λ r KL(r, p 0 ) + λ θ R θ (θ)where λ = {λ r , λ θ } are the regularization parameters and p 0 is a predefined reference distribution. As long as R θ is l.s.c., coercive and λ θ , λ r > 0, the resulting regularizer R E λ will also be l.s.c. and coercive. There are many possible choices for p 0 and R θ . For p 0 , this regularizes solutions of the gradient flow toward it; as such, in settings where there is some knowledge or preference about r at hand, we can set p 0 to reflect that. In our experiments, we utilize p 0 = N (0, M ) where M is a positive definite (p.d.) matrix. As for R θ , there are also many choices. In the context of neural networks, one natural choice is Tikhonov regularization 1 2 ∥ • ∥ 2 , resulting in weight decay (Hanson and Pratt, 1988) for gradient descent (Loshchilov and Hutter, 2019) which is a popular method for regularizing neural networks. In our experiments, we either use Tikhonov regularization or its simple variant ∥θ∥ 2 M := ⟨θ, M θ⟩.To solve the problem in Eq. ( 3), we construct a gradient flow that minimizes E λ . To this end, we endow the space M with a suitable notion of gradient∇ M E λ (θ, r) := (∇ θ E λ , ∇ r E λ )where ∇ θ and ∇ r denotes the Euclidean gradient and Wasserstein-2 gradient (Jordan et al., 1998) , respectively. The latter gradient is given by ∇r E λ (θ, r) := -∇ z • (r∇ z δ r E λ [θ, r]), where ∇ z • denotes the standard divergence operator and δ r denotes the first variation which is characterized in the following proposition.Proposition 4 (First Variation of E λ and RE λ ). Assume that E k θ (X|•) log q θ,r (X) p(X,y) < ∞ for all (θ, r) ∈ M; then the first variation of E λ is δ r E λ = δ r E +δ r R λ where δ r E[θ, r](z) = E k θ (X|z) log q θ,r (X) p(X,y) , and δ r R E λ [θ, r] = λ r log r/p 0 .The proof can be found in App. E.3. Thus, the (Euclidean-Wasserstein) gradient flow of E λ isEQUATIONWe now establish that the above gradient flow dynamic is contractive and that if a log-Sobolev inequality Eq. ( 7) holds, one can also establish exponential convergence. The log-Sobolev inequality is closely related to Polyak-Łojasiewicz inequality (or gradient dominance condition) and is commonly assumed in gradient-based systems to obtain convergence (for instance, see Kim et al. (2024) ). This is formally stated in the following proposition and proved in App. E.3. Proposition 5 (Contracting Gradient Dynamics). The free energy E λ along the flow Eq. ( 5) is non-increasing and satisfiesEQUATIONwhere∥∇ M E λ (θ, r)∥ 2 := ∥∇ θ E λ (θ, r)∥ 2 + ∥∇ z δ r E λ [θ, r]∥ 2 r .Moreover, if a log-Sobolev Inequality holds for a constant τ ∈ R >0 , i.e., for all (θ, r) ∈ M, we haveEQUATION)where E * λ := inf (θ,r)∈M E λ (θ, r); then we have exponential convergenceE λ (θ t , r t ) -E * λ ≤ exp(-t/τ )(E λ (θ 0 , r 0 ) -E * λ ).Typically direct simulation of the gradient flow Eq. ( 5) is intractable as the derivative of the first variation of R E λ involves ∇ z log r t ; instead, it is useful to identify the gradient flow with a McKean-Vlasov SDE, for which they share the same Fokker-Planck equation. The key distinction is that the SDE can be simulated without access to this quantity ∇ z log r t . This SDE, which we term the PVI flow, is given bydθ t = -∇ θ E λ (θ t , r t ) dt, dZ t = b(θ t , r t , Z t ) dt + 2λ r dW t , where r t = Law(Z t ), (8)where the drift is b(θ, r, •) := -∇ z δ r E[θ, r] + λ r ∇ z log p 0 (with the first variation given in Prop. 4) and W t is a d z -dimensional Wiener process. A connection between the Langevin diffusion, i.e., dZ t = ∇ z log p(Z t , y) dt + √ 2dW t , and PVI flow can be observed with the fixed kernel k θ (dx|z) = δ z (dx) and λ r = 0, namely, they both satisfy the same Fokker-Planck equation.Algorithm 1 Particle Variational Inference (PVI)Input: initialization (θ 0 , {Z 0 } M i=1 ); regularization parameters {λ θ , λ r }; step-sizes h θ and h r ; number of Monte Carlo samples L (for Eqs. ( 11) and ( 12)); and preconditioner Ψ = (Ψ θ , Ψ r ).for k = 1 to K do r M k-1 ← 1 M M m=1 δ Z k-1,m θ k ← θ k-1 -h θ Ψ θ ∇ θ E λ (θ k-1 , r M k-1 ) ▷ See Eq. (11) bk ← Z → -∇ z δ r E[θ k , r M k-1 ](Z) + λ r ∇ z log p 0 (Z) ▷ See Eq. (12) for m = 1 to M do Z k,m ← Z k-1,m + h r Ψ rb (Z k-1,m ) + √ λ r h r Ψ r η k,m ▷ η k,m ∼ N (0, I dz ) end for end for return (θ K , {Z K,m } M m=1 )To produce a practical algorithm, we are faced with several practical issues. The first issue we tackle is the computation of gradients of expectations for which using standard automatic differentiation is insufficient. The second problem is that these gradients are often ill-conditioned and have different scales in each dimension. This is tackled using preconditioning resulting in adaptive stepsize. Finally, to produce computationally feasible algorithms, we show how to discretize the PVI flow in both space and time. PVI is summarised in Algorithm 1.Computing the gradients. In the PVI flow, both the drift of the ODE and SDE include a gradient of an expectation with respect to parameters that define the distribution that is being integrated. Specifically, the terms that contain these gradients are ∇ θ E λ and ∇ z δ r E λ . Fortunately, these gradients can be rewritten as an expectation (as described in Prop. 6) for which the parameters being differentiated w.r.t. is only found in the integrand (see derivation in App. G.1). Proposition 6. If ϕ and k are differentiable, then we haveEQUATION)where ∇ θ ϕ ∈ R d θ ×dx denotes the Jacobian (∇ θ ϕ) ij = ∂ θi ϕ j (and similarly for ∇ z ϕ); scores are s θ,r (z, ϵ) := ∇ x log q θ,r (ϕ θ (z, ϵ)) (and similarly s p (z, ϵ)) ; and • denotes the usual matrix-vector multiplication in the sense ofM • v : (z, ϵ) → M (z, ϵ)v(z, ϵ).From Eqs. ( 9) and (10), we can produce Monte Carlo estimators for the gradients, i.e.,EQUATIONEQUATIONwhere {ϵ l } L l=1 i.i.d.∼ p k . This is an instance of a path-wise Monte-Carlo gradient estimator; a performant estimator that has been shown empirically to exhibit lower variance than other standard estimators (Kingma and Welling, 2014; Roeder et al., 2017; Mohamed et al., 2020) .Adaptive Stepsizes. One of the complexities of training neural networks is that their gradient is often poorly conditioned. As a result, for certain problems, the gradients computed from Eq. ( 9) and Eq. ( 10) can often produce unstable algorithms without careful tuning of the step sizes. When this occurs, we utilize preconditioners (Staib et al., 2019) to avoid this issue. Let Ψ θ : Θ → R d θ ×d θ and Ψ r : R dz → R dz×dz be the precondition for components θ and r respectively, then the resulting preconditioned gradient flow is given byEQUATION)If Ψ θ and Ψ r are positive definite, then E λ (θ t , r t ) remains non-increasing, i.e., Eq. ( 6) holds. As before, this Fokker-Planck equation is satisfied by the following Mckean-Vlasov SDE:EQUATION)EQUATION)where(∇ z • Ψ r ) i = dz j=1 ∂ zj [(Ψ r) ij ] and r t = Law(Z t ). The equivalence between Eq. ( 13) and Eqs. ( 14) and ( 15) is shown in App. G.2. A simple example for the preconditioner allows the θ t and Z t to have different time scales; ultimately, this results in different step sizes. Another more complex example of preconditioner Ψ θ is the RMSProp (Tieleman and Hinton, 2012) , and Ψ r we utilize a preconditioner inspired by RMSProp (see App. G.2). As with other related works (e.g., see Li et al. (2016) ), we found that the additional term ∇ z • Ψ r can be omitted in practice: it has little effect but incurs a large computational cost. Discretization in both space and time. To obtain an actionable algorithm, we need to discretize the PVI flow in both space and time. For the space discretization, we propose to use a particle approximation for r t , i.e., for a set of particles {Z t,m } M m=1 with each satisfying Law(Z t,m ) = r t , we utilize the approximation r M t := 1 M M m=1 δ Zt,m which converges almost surely to r t in the weak topology as M → ∞ by the strong law of large numbers and a countable determining class argument (e.g., see Schmon et al. (2020, Theorem 1.1) ). This approximation is key to making the intractable tractable, e.g., Eq. ( 1) is approximated by q θ,r M t = 1 M M m=1 k θ (x|Z t,m ). One obtains a particle approximation to the PVI flow from the following ODE-SDE:dθ M t = -∇ θ E λ (θ M t , r M t ) dt, ∀m ∈ [M ] : dZ M t,m = b(θ M t , r M t , Z M t,m ) dt + 2λ r dW t,m, where [M ] := {1, . . . , M }. As for the time discretization, we employ Euler-Maruyama discretization with step-size h which (using an appropriately defined preconditioner) can be decoupled into different stepsizes for θ t and Z t denoted by h θ and h r respectively.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
theoretical analysis: 
 We are interested in the behaviour of the PVI flow (8). However, a key issue in its study is that the drift in PVI flow might lack the necessary continuity properties to analyze using the existing theory. In this section, we instead analyze the related gradient flow of the more regular functionalEQUATIONfor γ > 0. A similar modified flow was also explored in Crucinio et al. (2024) for similar reasons; they found empirically that, at least when using a tamed Euler scheme, setting γ = 0 did not cause problems in practice. Similarly, our experimental results for PVI found γ = 0 did not have issues. To provide an additional measure of confidence in the reasonableness of this regularization and of the use of this functional as a proxy for E λ , we establish that the minima of E γ λ converge to those of E λ in the γ → 0 limit. Proposition 7 (Γ-convergence and convergence of minima). Under the same assumptions as Prop. 3, we have that E γ λ Γ-converges to E λ as γ → 0 (in the sense of Def. B.1). Moreover, we have as an immediate corollary thatinf (θ,r)∈M E λ (θ, r) = lim γ→0 inf (θ,r)∈M E γ λ (θ, r).The proof uses techniques from Γ-convergence theory (introduced by De Giorgi, see e.g. De Giorgi and Franzoni (1975) ; see Dal Maso (2012) for a good modern introduction) and can be found in App. F.1. The gradient flow of E γ λ , which we term γ-PVI, is given byEQUATIONwherer γ t = Law(Z γ t ) and b γ (θ, r, •) = -∇ z δ r E γ [θ, r] + λ r ∇ z log p 0 .The derivation follows similarly to that in Section 3.2, and is omitted for brevity. The use of γ > 0 is crucial for establishing key regularity conditions in our analysis. We proceed by stating our assumptions. Assumption 1 (Regularity of the target p, reference distribution p 0 , and R θ ). We assume that log p(y) is bounded; and p, p 0 and R θ have Lipschitz gradients with constants K p , K p0 , K R θ respectively: there exists some B p ∈ R >0 such that log p(y) ≤ B p ; and for any given y there exists a K p ∈ R >0 such that ∥∇ x log p(x, y) -∇ x log p(x ′ , y)∥ ≤ K p ∥x -x ′ ∥ for all x, x ′ ∈ R dx (similarly for p 0 and R θ ). Assumption 2 (Regularity of k). We assume that the kernel k and its gradient is bounded and has K k -Lipschitz gradient; i.e., there exist constantsB k , K k ∈ R >0 such that |k θ (x|z)| + ∥∇ (θ,x,z) k θ (x|z)∥ ≤ B k , and ∥∇ x k θ (x|z) -∇ x k θ ′ (x ′ |z ′ )∥ ≤ K k (∥(θ, x, z) -(θ ′ , x ′ , z ′ )∥) hold for all θ, θ ′ ∈ Θ, z, z ′ ∈ R dz ,and x, x ′ ∈ R dx , Assumption 3 (Regularity of ϕ and p k .). We assume that ϕ has Lipschitz gradient and bounded gradient. In order words, there isa ϕ ∈ R ≥0 , b ϕ ∈ R >0 such that ϕ satisfies ∥∇ (θ,z) ϕ(z, ϵ) - ∇ (θ,z) ϕ θ ′ (z ′ , ϵ)∥ F ≤ (a ϕ ∥ϵ∥ + b ϕ )(∥(θ, z) -(θ ′ , z ′ )∥) and ∥∇ (θ,z) ϕ θ (z, ϵ)∥ F ≤ (a ϕ ∥ϵ∥ + b ϕ ) for all (θ, z), (θ ′ , z ′ ) ∈ Θ × R dz , ϵ ∈ R dx ,where ∥ • ∥ F denotes the Frobenius norm. We also assume that p k has finite second moments.Assumptions 2 and 3 are intimately connected; under some regularity conditions, one may imply the other but we shall abstain from this digression for the sake of clarity. Assumptions 2 and 3 are quite mild and hold for popular kernels such as k θ (x|z) = N (x; µ θ (z), Σ) under some regularity assumptions on µ θ and Σ (which we show in App. C). These assumptions are key to establishing that the drift in Eq. ( 17) is Lipschitz continuous (see Prop. 11 in the App.), from which, we establish the existence and uniqueness of the solutions of Eq. ( 17). Proposition 8 (Existence and Uniqueness). Under Assumptions 1 to 3, if γ > 0 andE p k (ϵ) ∥s γ θ,r (z, ϵ) -s p (z, ϵ)∥ is bounded (with s γ θ,r : (z, ϵ) → ∇ x log(q θ,r • ϕ θ (z, ϵ) + γ)); then, given (θ 0 , r 0 ) ∈ M, the solutions to Eq. (17) exists and is unique.The proof can be found in App. F.2. Under the same assumptions, we can establish an asymptotic propagation of chaos result that justifies the usage of a particle approximation in place of r γ t in Eq. ( 17). Proposition 9 (Propagation of chaos). Under the same assumptions as Prop. 8; we have for any fixed T :lim M →∞ E sup t∈[0,T ] θ γ t -θ γ,M t 2 + W 2 2 (r γ t ) ⊗M , q γ,M t = 0, where (r γ t ) ⊗M = M i=1 (r γ t ); q γ,M t = Law({Z γ,M t,m } M m=1 ); θ γ,M t and Z γ,M t,m are solutions to dθ γ,M t = -∇ θ E γ λ (θ γ,M t , r γ,M t ) dt, where r γ,M t = 1 M M m=1 δ Z γ,M t,m ∀m ∈ [M ] : dZ γ,M t,m = b γ (θ γ,M t , r γ,M t , Z γ,M t,m ) dt + 2λ r dW t,m .The proof can be found in App. F.3. Having established the existence and uniqueness of the γ-PVI flow, as well as an asymptotic justification for using particles, we now provide a numerical evaluation to demonstrate the efficacy of our proposal. 
experiments: 
 In this section, we compare PVI against other semi-implicit VI methods. As described in the App. A, these include unbiased semi-implicit variational inference (UVI) of Titsias and Ruiz (2019), semiimplicit variational inference (SVI) of Yin and Zhou (2018) , and the score matching approach (SM)of Yu and Zhang (2023) . Through experiments, we show the benefits of optimizing the mixing distribution; we compare the effectiveness of PVI against other SIVI methods on a density estimation problem on toy examples; and, we compare against other SIVI methods on posterior estimation tasks for (Bayesian) logistic regression and (Bayesian) neural network. The details for reproducing experiments as well as computation information can be found in App. H. The code is available at https://github.com/jenninglim/pvi. The plot shows the density q θ,r from PVI and PVIZero as well as the KDE plot of r from PVI described by 100 particles.From Prop. 1, it can be said that ultimately current SIVI methods utilize (directly or indirectly) a fixed mixing distribution whilst PVI does not. We are interested in establishing whether there is any benefit to optimizing the mixing distribution. Intuitively, the mixing distribution can be utilized to express complex properties, such as multimodality, which the neural network kernel k θ can then exploit. If the mixing distribution is fixed, this means that the neural network must learn to express these complex properties directly-which can be difficult (Salmona et al., 2022) . This intuition turns out to hold, but for the kernel to exploit an expressive mixing distribution, it must be designed well. To illustrate this, consider the distributions 1 2 N (µ, I) + 1 2 N (-µ, I) for µ = {1, 2, 4} and following kernels: the "Constant" kernel N (z, I 2 ); "Push" kernelN (f θ (z), σ 2 θ I 2 ); "Skip" kernel N (z + f θ (z), σ 2 θ I 2 ); and "LSkip" kernel N (W z + f θ (z), σ 2 θ I 2 ) where W ∈ R 2×2 ;. We compare the results from PVI and PVIZero (PVI with h r = 0 to result in a fixed r ≈ N (0, I 2 )) to emulate PVI with a fixed mixing distribution. As µ gets larger, the complexity of the kernel (or the mixing distribution) must grow to express this (e.g., see (Salmona et al., 2022 , Corollary 2)).Fig. 1 shows the resulting densities and the learnt mixing distribution of PVI and PVIZero for different kernels and various µ. For the constant kernel, PVI can solve this problem by learning a complex mixing distribution to express the multimodality. However, for the push kernel, it can be seen that as µ gets larger PVI and PVIZero suffer from mode collapse which we suspect is due to the mode-seeking behaviour of using reverse KL and why prior SIVI methods utilized annealing methods (see Yu and Zhang (2023, Section 4.1) ). As a remedy, we utilize a Skip kernel which can be seen to improve both PVI and PVIZero. In particular, both PVI and PVIZero were able to successfully express the bimodality in µ = 2; however, PVIZero falls short when µ = 4 while PVI can express the multimodality by learning a bimodal mixing distribution. Since Skip requires d z = d x , we show that LSkip (which removes the requirement) exhibits a similar behaviour to Skip.We follow prior works (e.g., Yin and Zhou (2018) ) and consider three toy examples whose densities are shown in Fig. 2 (they are given explicitly in App. H.2). In this setting, we use the kernelk θ (x|z) = N (x; z + f θ (z), σ 2 θ I) with d z = d x = 2where f θ (z) is a neural network whose architecture can be found in App. H.2. As a qualitative measure of performance, Fig. 2 shows the resulting approximating distribution of PVI which can be seen to be a close match to the desired distribution. To compare methods quantitatively, we report the (sliced) Wasserstein distance Figure 2 : Contour plots of the densities q θ,r (in blue) against the true densities (in black) for various toy density estimation problems. We also plot the absolute difference in the density of q θ,r and the true density, i.e., |q θ,r -p|.(computed by POT (Flamary et al., 2021) ) and the rejection power of a state-of-the-art two-sample kernel test (Biggs et al., 2023) between the approximating and true distribution in Table 1 . The results reported are the average and standard deviation (from ten independent trials of the respective SIVI algorithms). In each trial, the rejection rate p is computed from 100 tests and the sliced Wasserstein distance is computed from 10000 samples with 100 projections. If the variational approximation matches the distribution, the rejection rate will be at the nominal level of 0.05. It can be seen that PVI consistently performs better than SIVI across all problems. PVI can achieve a rejection rate near nominal levels across all problems whilst other algorithms can achieve good performances on one but not the other. The details regarding how the Wasserstein distance is calculated and the hyperparameters used can be found in App. showing the standard deviation estimated from 10 independent runs. We indicate in bold when the rejection rate minus the standard deviation is lower than the nominal level 0.05, and for the algorithm that achieves the lowest Wasserstein score.As with others (Yin and Zhou, 2018) , we consider a Bayesian logistic regression problem on the waveform dataset (Breiman and Stone, 1984) . The model is expressed as y | x, w ∼ Bernoulli(Sigmoid(⟨x, w⟩)) with prior x ∼ N (0, 0.01 -1 × I 22 ) where (y, w) ∈ {0, 1} × R 21 is the response and covariates, and w := [1, w] is the covariates with appended one for the intercept. The "ground truth" is composed of posterior samples generated from running Markov chain Monte Carlo (MCMC) samples in Yin and Zhou (2018) . We use the kernelk θ (x|z) = N (W z + f θ (z), exp( 1 2 [M θ + M ⊤ θ ]))where exp denotes the matrix exponential which ensures positive definiteness. In Fig. 3 , we visually compare certain statistics of the distribution obtained from MCMC and distribution obtained from SIVI methods. Fig. 3a shows the pair-wise marginal posterior distributions for three weights x 1 , x 2 , x 3 chosen at random from MCMC and SIVI approximations; and in Fig. 3b we compare the correlation coefficients obtained from MCMC against SIVI methods. It can be seen that PVI obtains an approximation close to MCMC with most other SIVI methods obtaining similar performance levels (with the exception of SM). See App. H.3 for all the implementation details.Following prior works (e.g., Yu and Zhang ( 2023)), we compare our methods with other baselines on sampling the posterior of the Bayesian neural network for regression problems on a range of real-world datasets. We utilize the LSkip kernel (Gerritsma et al., 2013) 0.13 0.02 0.18 0.02 0.17 0.02 0.98 0.16 Table 2 : Root mean square error (lower is better) for Bayesian neural networks on the test set for various datasets. Here, we write the results in the form µ σ where µ is the average RMS and σ is its standard error computed over 10 independent trials. We indicate in bold the lowest score.k θ (x|z) = N (x; W z + f θ (z), σ 2 θ (z)I dx ). In 
conclusion, limitations, and future work: 
 In this work, we frame SIVI as a minimization problem of E λ , and then, as a solution, we study its gradient flow. Through discretization, we propose a novel algorithm called Particle Variational Inference (PVI). Our experiments found that PVI can outperform current SIVI methods. At a marginal increase in computation cost (see App. H) compared with prior methods, PVI can consistently perform better (or at least comparably in the worst cases considered) which we attribute to not imposing a particular form on the mixing distribution. This is a key advantage of PVI compared to prior methods: by not relying upon push-forward mixing distributions and instead using particles, the mixing distribution can express arbitrary distributions when the number of particles is sufficiently large. Furthermore, it is not necessary to tune the family of mixing distributions to obtain good results in particular problems. Theoretically, we study a related gradient flow of E γ λ and establish desirable properties such as the existence and uniqueness of solutions and propagation of chaos results.The main limitation of our work is that the theoretical results only apply to the case where γ > 0; yet, our experiments were performed with γ = 0 as this is when E λ corresponds to the (regularized) evidence lower bound. In future work, one can address these limitations by reducing this gap. Furthermore, we found that certain kernels were more amenable than others when exploiting an expressive mixing distribution (e.g., the skip kernel). The question of designing these kernels for PVI (or SIVI more generally) is important for future work. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: We do not release data or models that have a high risk of misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: Yes, we reference the dataset used.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: In this section, we outline four areas of related work: semi-implicit variational inference; Euclidean-Wasserstein gradient flows and Wasserstein-gradient flows in VI; mixture models in VI; and finally, the link between SIVI and solving Fredholm equations of the first kind.At the time of writing, there are three algorithms for SIVI proposed: SVI Yin and Zhou (2018) , UVI Titsias and Ruiz (2019), SM Yu and Zhang (2023) . Concurrently, Cheng et al. (2024) extended the SM variant by solving the inner minimax objective and simplified the optimization problem. Each had their parameterization of SID (as discussed in Section 2), and their proposed optimization method. SVI relies on optimizing a bound of the ELBO which is asymptotically tight. UVI, like our approach, optimizes the ELBO by using gradients-based approaches. However, one of its terms is the score ∇ x log q θ,r (x) which is intractable. The authors proposed using expensive MCMC chains to estimate it; in contrast to PVI, this term is readily available to us. For SM, they propose to optimize the Fisher divergence, however, to deal with the intractabilities the resulting objective is a minimax optimization problem which is difficult to optimize compared to standard minimization problems.PVI utilizes the Euclidean-Wasserstein geometry. This geometry and associated gradient flows are initially explored in the context of (marginal) maximum likelihood estimation by Kuntz et al. (2023) and their convergence properties are investigated by Caprio et al. (2024) . In Lim et al. (2024) , the authors investigated accelerated gradient variants of the aforementioned gradient flow in Euclidean-Wasserstein geometry. The Wasserstein geometry particularly for gradient flows on probability space has received much attention with many works exploring different functionals (for examples, see Arbel et al. (2019) ; Korba et al. (2021) ; Li et al. (2023) ). In the context of variational inference Lambert et al. (2022) analyzed VI as a Bures-Wasserstein Gradient flow on the space of Gaussian measures.PVI is reminiscent of mixture distributions which is a consequence of the particle discretization. Mixture models have been studied in prior works as variational distributions (Graves, 2016; Morningstar et al., 2021) . In Graves (2016) , the authors extended the parameterization trick to mixture distributions; and Morningstar et al. (2021) proposed to utilize mixture models as variational distributions in the framework of Kucukelbir et al. (2017) . Although similar, the mixing distribution assists the kernel in expressing complex properties of the true distribution at hand (see Section 5.1) which is an interpretation that mixture distribution lacks.There is an obvious similarity between SIVI and solving Fredholm equations of the first kind. There is considerable literature on solving such problems; see Crucinio et al. (2024) , which is closest in spirit to the approach of the present paper, and references therein. In fact, writing p(•|y) = k(•|z, θ)r(z)dz. with k(•|z, θ) ≡ k θ (•|z) makes the connection more explicit: essentially, one seeks to solve a nonstandard Fredholm equation, with the LHS known only up to a normalizing constant, constraining the solution to be in P(Z) × {δ θ : θ ∈ Θ}. While Crucinio et al. (2024) develop and analyse a simple Wasserstein gradient flow to address a regularised Fredholm equation, neither the method nor analysis can be applied to the SIVI problem because of this non-trivial constraint.The following is one of many essentially equivalent definitions of Γ-convergence (see Dal Maso (2012); Braides (2002) for comprehensive summaries of Γ-convergence). We take as definition the following (see Dal Maso (2012, Proposition 8.1), Braides (2002, Definition 1.5)): Definition B.1 (Γ-convergence). Assume that M is a topological space that satisfies the first axiom of countability. Then a sequence F γ : M → R is said to Γ-converge to F if:• (lim-inf inequality) for every sequence (θ γ , r γ ) ∈ M converging to (θ, r) ∈ M, we havelim inf γ→0 F γ (θ γ , r γ ) ≥ F(θ, r).• (lim-sup inequality) for any (θ, r) ∈ M, there exists a sequence (θ γ , r γ ) ∈ M, known as a recovery sequence, converging to (θ, r) which satisfieslim sup γ→0 F γ (θ γ , r γ ) ≤ F(θ, r).Γ-convergence corresponds, roughly speaking, to the convergence of the lower semicontinuous envelope of a sequence of functionals and, under mild further regularity conditions such as equicoercivity, is sufficient to ensure the convergence of the sets of minimisers of those functionals to the set of minimisers of the limit functional.We shall show that the Gaussian kernel k θ (x|z) = N (x; µ θ (z), Σ), i.e.,k θ (x|z) = (2π) -dx/2 det(Σ) -0.5 exp - 1 2 (x -µ θ (z)) T Σ -1 (x -µ θ (z)) ,where µ θ : R dz → R dx ; and Σ ∈ R dx×dx and is positive definite. In this section, we show that Assumptions 2 and 3 are implied by Assumptions 4 and 5. Assumption 4. µ θ is bounded and Σ is positive definite: there exists B µ ∈ R >0 such that the following holds for all (θ, z)∈ Θ × R dz : ∥∇ (θ,z) µ θ (z)∥ F ≤ B µ ,and for any x ∈ R dx \ 0, x T Σx > 0. Assumption 5. µ θ is Lipschitz and has Lipschitz gradient, i.e., there exist constants K µ ∈ R >0 such that for all (θ, z), (θ ′ , z ′ ) ∈ Θ × R dz the following hold:∥µ θ (z) -µ θ ′ (z ′ )∥ ≤ K µ ∥(θ, z) -(θ ′ , z ′ )∥, ∥∇ (θ,z) µ θ (z) -∇ (θ,z) µ θ ′ (z ′ )∥ F ≤ K µ ∥(θ, z) -(θ ′ , z ′ )∥.C.1 k θ satisfies Assumption 2In this section, we show that k θ satisfies Assumption 2. We first show the boundedness property then the Lipschitz property.Boundedness. First, we shall show that k θ is bounded. Clearly, we have k θ (x|z) ∈ 0, (2π) -dx/2 det(Σ) -0.5 hence |k θ | is bounded as a consequence of Assumption 4. Now to show that the gradient is bounded ||∇ (θ,x,z) k θ (x|z)||, we have the following∇ x k θ (x|z) = -k θ (x|z)Σ -1 (x -µ θ (z)), ∇ z k θ (x|z) = ∇ z µ θ (z)∇ µ N (x; µ, σ 2 I dx ) | µ θ (z) , ∇ θ k θ (x|z) = ∇ θ µ θ (z)∇ µ N (x; µ, σ 2 I dx ) µ θ (z) .Hence, we have ∥∇ (x,µ,σ) N (x; µ θ (z), σ 2 I dx )∥ < ∞, from Assumption 4 and using the fact the gradient of a Gaussian density of given covariance w.r.t. µ is uniformly bounded. Thus, we have shown that k θ satisfies the boundedness property in Assumption 2.Lipschitz. For k θ , one choice of coupling function and noise distribution is ϕ θ (z, ϵ) = Σ 1 2 ϵ + µ θ (z) and p k = N (0, I dx ) where Σ 1 2 be the unique symmetric and positive definite matrix with (Σ 1 2 ) 2 = Σ (Horn and Johnson, 2012, Theorem 7.2.6); and the inverse map is ϕ-1 θ (z, x) = Σ -1 2 (x -µ θ (z)). Thus, from the change-of-variables formula, we have∇ x k θ (x|z) = ∇ x [p k (ϕ -1 θ (z, x))det(∇ x ϕ -1 θ (z, x))] = det(∇ x ϕ -1 θ (z, x))∇ x [p k ϕ -1 θ (z, x) ] = det(Σ -1/2 )Σ -1/2 ∇ x p k (ϕ -1 θ (z, x)) = Σ-1 2 ∇ x p k (ϕ -1 θ (z, x)) where Σ-1 2 := det(Σ -1/2 )Σ -1/2. Thus, we have∥∇ x k θ (x|z) -∇ x k θ ′ (x ′ |z ′ )∥ ≤ ∥ Σ-1 2 ∇ x p k (ϕ -1 θ (z, x)) -Σ-1 2 ∇ x p k (ϕ -1 θ ′ (z ′ , x ′ ))∥ ≤ ∥ Σ-1 2 ∥ F ∥∇ x p k (ϕ -1 θ (z, x)) -∇ x p k (ϕ -1 θ ′ (z ′ , x ′ ))∥ ≤ C∥ϕ -1 θ (z, x) -ϕ -1 θ ′ (z ′ , x ′ )∥,where C is a constant and we use the following facts: ∥ Σ-1 2 ∥ F ≤ |det(Σ -1/2 )|∥Σ -1/2 ∥ F < ∞ following from the fact Σ -1/2 is positive definite; p k is a standard Gaussian density function with Lipschitz gradients; and that the inverse map ϕ -1 is Lipschitz from Assumptions 4 and 5:∥ϕ -1 θ (z, x) -ϕ -1 θ ′ (z ′ , x ′ )∥ ≤ ∥Σ -1 2 ∥ F ∥(x, µ θ (z)) -(x ′ , µ θ ′ (z ′ ))∥ ≤ C ′ ∥(x, θ, z) -(x ′ , θ ′ , z ′ )∥.Hence, we have shown that k θ satisfies the Lipschitz property of Assumption 2, and so Assumption 2 holds for k θ .C.2 k θ satisfies Assumption 3One can compute the gradient as∇ (θ,z) ϕ θ (z, ϵ) = ∇ (θ,z) µ θ (z),and hence ∥∇ (θ,z) ϕ θ (z, ϵ)∥ F is bounded from Assumption 4. The Lipschitz gradient property is immediate from Assumption 5.p k has finite second moments since it is a Gaussian.Proof of Prop. 1. We start by showing Q YuZ = Q TR . To this end, we begin by showing the inclusionQ YuZ ⊆ Q TR , i.e., Q(K F ;ϕ,p k , R G,pr ) ⊆ Q(K F •G;ϕ,p k , {p r }). Let q ∈ Q(K F ;ϕ,p k , R G,pr), then there is some f ∈ F and g ∈ G such that q = q k f ;ϕ,p k ,g # pr . From straight-forward computation, we haveq k f ;ϕ,p k ,g # pr = E z∼g # pr [k f ;ϕ,p k (•|z)] (a) = E z∼pr [k f ;ϕ,p k (•|g(z))] ∈ Q(K F •G;ϕ,p k , {p r }),where (a) follows the law of the unconscious statistician, and the last element-of follows from the fact that kf ;ϕ,p k (•|g(ϵ))) = ϕ(f • g(ϵ), •) # p k ∈ K F •G;ϕ,p k .We can follow the argument above in reverse to obtain the reverse inclusion. Hence, we have obtained as desired.That Q YuZ = Q YiZ , follows in a similar manner, which we shall outline for completeness: let q ∈ Q(K F ;ϕ,p k , R G,pr ), thenq = q k f ;ϕ,p k ,g # pr = E z∼g # pr [k f ;ϕ,p k (•|z)] = E z∼f •g # pr [k ϕ,p k (•|z)] ∈ Q YiZ .One can conclude by applying the same logic in the reverse direction.E Proofs in Section 3E.1 Proof of Prop. 2Proof of Prop. 2. (E is lower bounded). Clearly, we haveE(θ, r) = KL(q θ,r , p(•|y)) -log p(y) ≥ -log p(y),Hence, we have E(θ, r) ∈ [-log p(y), ∞) which is lower bounded by our assumption.(E is lower semi-continuous). Let (θ n , r n ) n∈N be such that lim n→∞ r n = r and lim n→∞ θ n = θ.We can split the domain of integration, and write E equivalently asEQUATION+ 1 [0,1) p(x, y) q θ,r (x) log q θ,r (x) p(x, y)EQUATIONWe shall focus on the RHS of (18).Note that we have the following bound -1 [1,∞) p(x, y) q θn,rn (x) log q θn,rn (x) p(x, y) q θn,rn (x)≤ max 0, log p(x, y) q θn,rn (x) q θn,rn (x) ≤ max{0, C},where C is some constant. The last inequality follows from the fact that the evidence is bounded from above and the kernel is bounded. We can apply Reverse Fatou's Lemma to obtainlim sup n→∞ -1 [1,∞)p(x, y) q θn,rn (x) log q θn,rn (x) p(x, y) q θn,rn (x) dx≤ lim sup n→∞ -1 [1,∞)p(x, y) q θn,rn (x) log q θn,rn (x) p(x, y) q θn,rn (x) dx.Since we have the following relationshipslim sup n→∞ 1 [1,∞) p(x, y) q θn,rn (x) ≤ 1 [1,∞) p(x, y) q θ,r (x) , lim n→∞ -log p(x, y) q θn,rn (x) = -log p(x, y) q θ,r (x) , lim n→∞ q θn,rn = q θ,r pointwise,where the first line is from u.s.c. of 1 [1,∞) ; the second line from the continuity of log; the final line follows from the bounded kernel k assumption and dominated convergence theorem.Thus, we have thatlim sup n→∞ -1 [1,∞)p(x, y) q θn,rn (x) log q θn,rn (x) p(x, y) q θn,rn (x) dx≤ -1 [1,∞)p(x, y) q θ,r (x) log q θ,r (x) p(x, y) q θ,r (x) dx,Using the fact that lim sup n→∞ -x n = -lim inf n→∞ x n , we have shown that-lim inf n→∞ 1 [1,∞)p(x, y) q θn,rn (x) log q θn,rn (x) p(x, y) q θn,rn (x) dxEQUATIONSimilarly, for the RHS of (19), using Fatou's Lemma (with varying measure and the set-wise convergence of q θn,rn ) and using the l.s.c. of 1 [0,1) , we obtain that lim inf n→∞ 1 [0,1) p(x, y) q θn,rn (x) log q θn,rn (x) p(x, y) q θn,rn (x) dxEQUATIONHence, combining the bounds ( 20) and ( 21), we have that shown thatlim inf n→∞ E(θ n , r n ) = lim inf n→∞ log q θn,rn (x) p(x, y) q θn,rn (x) dx ≥ log q θ,r (x) p(x, y) q θ,r (x) dx ≥ E(θ, r).In other words, E is lower semi-continuous.(Non-Coercivity) To show non-coercivity, we will show that there exists some level set {(θ, r) : E(θ, r) ≤ β} that is not compact. We do this by finding a sequence contained in the level set that does not contain a (weakly) converging subsequence.Consider the sequence Π := (θ n , r n ) n∈N where θ n = θ 0 ; ∥θ 0 ∥ < ∞; r n = δ n ; k θ (x|z) = N (x; θ, I dx ); and p(x|y) = N (x; 0, I dx ). Clearly, we have q θ,r (x) = N (x; θ, I dx ) and so KL(q θ,r , p(•|y)) = 1 2 ∥θ∥ 2 . Hence, there is a β < ∞ such thatE(θ n , r n ) = KL(q θn,rn , p(•|y)) -log p(y) ≤ 1 2 ∥θ 0 ∥ 2 -log p(y) ≤ β.Thus, we have shown that Π ⊂ {(θ, r) : E(θ, r) ≤ β}. However, since the support of the elements of {r n ∈ P(R dz )} n eventually lies outside a ball of radius R for any R < ∞ and hence of any compact set, Π is not tight. Hence, Prokhorov's theorem (Shiryaev, 1996, p. 318 ) tells us that, as Π is not tight, it is not relatively compact. We conclude that, as the level set is not relatively compact, the functional is not-coercive.Proof of Prop. 3. (Coercivity) Consider the level set {(θ, r) : E λ (θ, r) ≤ β}, which is contained in a relatively compact set. To see this, first note that {(θ, r) : E λ (θ, r) ≤ β} ⊆ {(θ, r) : -log p(y) + R λ (θ, r) ≤ β} ⊆ {(θ, r) : R λ (θ, r) ≤ β + log p(y)} By coercivity of R λ , i.e., the above level set is relatively compact hence E λ is coercive.(Lower semi-continuity) Lower semi-continuity (l.s.c.) follows immediately from the l.s.c. of E and R λ .(Existence of a minimizer) The existence of a minimizer follows from Dal Maso (2012, Theorem 1.15) utilizing coercivity and l.s.c. of E λ .Recall from Santambrogio (2015, Definition 7.12), Definition E.1 (First Variation). If p is regular for F , the first variation of F : P(R dz ) → R, if it exists, is the element that satisfieslim ϵ→0 F (p + ϵχ) -F (p) ϵ = δ r F [r](z)χ(dz),for any perturbation χ = p -p with p ∈ P(R dz ) ∩ L ∞ c (R dz ) (see Santambrogio (2015, Notation) ). One can decompose the first variation of E γ λ as:δ r E γ λ [θ, r] = δ r E γ [θ, r] + δ r R E λ [θ, r]. where E γ : (θ, r) → log q θ,r (x)+γ p(x,y) q θ,r (dx). Since δ r R E λ [θ, r] = λ r δ r KL(r|p 0 ), its first variation follows immediately from standard calculations (Ambrosio et al., 2005; Santambrogio, 2015) . As for δ r E γ , we have the following proposition: Proposition 10 (First Variation of E γ ). Assume that for all (θ, r, z) ∈ M × R dz ,E k θ (X|z) log q θ,r (X) + γ p(X, y) < ∞,then we obtainδ r E γ [θ, r](z) = E k θ (X|z) log q θ,r (X) + γ p(X, y) + q θ,r (X) q θ,r (X) + γ .Proof. Since q θ,r+ϵχ = k θ (•|z)(r + ϵχ)(z) dz = q θ,r + ϵq θ,χ , we haveE γ (θ, r + ϵχ) = X q θ,r+ϵχ (x) log q θ,r+ϵχ (x) + γ p(y, x) dx = X [q θ,r + ϵq θ,χ ](x) log([q θ,r + ϵq θ,χ ](x) + γ) dx - X [q θ,r + ϵq θ,χ ](x) log p(y, x) dx.Applying Taylor's expansion, we obtain (x + ϵy) log(x + ϵy + γ) = x log(x + γ) + ϵy log(x + γ) + x x+γ + o(ϵ), we obtainE γ (θ, r + ϵχ) = X q θ,r (x) log q θ,r (x) + γ p(y, x) dx + ϵ X q θ,χ (x) log q θ,r (x) + γ p(y, x) + q θ,r (x) q θ,r (x) + γ dx + o(ϵ).Hence, we obtainlim ϵ→ E γ (θ, r + ϵχ) -E γ (θ, r) ϵ = X q θ,χ (x) log q θ,r (x) + γ p(y, x) + q θ,r (x) q θ,r (x) + γ dx = X Z k θ (x|z)χ(dz) log q θ,r (x) + γ p(y, x) + q θ,r (x) q θ,r (x) + γ dx (a) = Z X k θ (x|z) log q θ,r (x) + γ p(y, x) + q θ,r (x) q θ,r (x) + γ dx χ(z) dz.One can then identify the desired result. In (a), we appeal to Fubini's theorem for the interchange of integrals whose conditionsEQUATIONare satisfied by our assumptions. This can be seen from LHS Eq. ( 22)≤ Z E k θ (X|z) log q θ,r (X) + γ p(X, y) + q θ,r (X) q θ,r (X) + γ |χ(z)|dz ≤ 0,where we use our assumption and the fact that χ is absolutely integrable.Proof of Prop. 5. The result can be obtained from direct computation. We begind dt E λ (θ t , r t ) = ∇ θ E λ (θ t , r t ), θt + δ r E λ [θ t , r t ] ∂ t r t dzThe second term can be simplifiedδ r E λ [θ t , r t ] ∂ t r t dz = δ r E λ [θ t , r t ](z)∇ z • (r t (z)∇δ r E λ [θ t , r t ](z)) dz = -r t ∥∇ z δ r E λ [θ t , r t ](z)∥ 2 dzwhere the last inequality follows from integration by parts. Hence, the claim holds. If the log-Sobolev inequality holds, then we haved dt [E λ (θ t , r t ) -E * λ ] = -∥∇ M E λ [θ t , r t ]∥ ≤ - 1 τ [E λ (θ t , r t ) -E * λ ]. From Grönwalls inequality, we obtain the desired result.F.1 Proof of Prop. 7Proof. We first begin by proving Γ-convergence directly via its definition, i.e., demonstrating that the liminf inequality holds and establishing the existence of a recovery sequence. The latter follows from pointwise convergence:lim γ→0 E γ λ (θ, r) = E λ (θ, r), upon taking (θ γ , r γ ) = (θ, r) for all γ.The liminf inequality can be seen to follow similarly from the l.s.c. argument in App. E.1.To arrive at the convergence of minima, we invoke Dal Maso (2012, Theorem 7.8) by using the fact that E γ λ is equi-coercive in the sense of Dal Maso (2012, Definition 7.6). To see that E γ λ is equi-coercive, note that we have E γ λ ≥ E λ and E λ is l.s.c. (from Prop. 3), then applying Dal Maso (2012, Proposition 7.7).Proof of Prop. 8. We can equivalently write the γ-PVI flow in Eq. ( 17) as followsEQUATIONwhere σ = 0 0 0 √ 2λ r I dz , and bγ : Rd θ × P(Z) × R dz → R d θ +dz : (θ, r, Z) → -∇ θ E γ λ (θ, r) b γ (θ, r, Z) .In App. F.4, we show that under our assumptions the drift bγ is Lipschitz. And under Lipschitz regularity conditions, the proof follows similarly to Lim et al. (2024) which we shall outline for completeness.We begin endowing the space Θ × P(R dz ) with the metricd((θ, r), (θ ′ , r ′ )) = ∥θ -θ ′ ∥ 2 + W 2 2 (q, q ′ ).Let Υ ∈ C([0, T ], Θ × P(R dz )) and denote Υ t = (ϑ Υ t , ν Υ t ) for it's respective components. Consider the process that substitutes Υ into (23), in place of the Law(Z t ) and θ t ,d(θ Υ t , Z Υ t ) = bγ (ϑ Υ t , ν Υ t , Z Υ t ) dt + σ dW t .whose existence and uniqueness of strong solutions are given by Carmona (2016) [Thereom 1.2].F T : C([0, T ], Θ × P(R dz )) → C([0, T ], Θ × P(R dz )) : Υ → (t → (θ Υ t , Law(Z Υ t )). Let (θ t , Z t ) denote a process that is a solution to (23) then the function t → (θ t , Law(Z t )) is a fixed point of the operator F T . The converse also holds. Thus, it is sufficient to establish the existence and uniqueness of the fixed point of the operator F T . For Υ = (ϑ, ν) andΥ ′ = (ϑ ′ , ν ′ ) ∥θ Υ t -θ Υ ′ t ∥ 2 + E[∥Z Υ t -Z Υ ′ t ∥] 2 = E t 0 bγ (ϑ s , ν s , Z Υ s ) -bγ (ϑ ′ s , ν ′ s , Z Υ ′ s ) ds 2 ≤ tC t 0 E∥Z Υ s -Z Υ ′ s ∥ 2 + ∥ϑ s -ϑ ′ s ∥ 2 + W 2 1 (ν s , ν ′ s ) ds ≤ C(t) t 0 [W 2 2 (ν s , ν ′ s ) + ∥ϑ s -ϑ ′ s ∥ 2 ] ds,where we apply Jensen's inequality; C r -inequality; Lipschitz drift of bγ ; and Grönwall's inequality. The constant C := 3K 2 b and C(t) := tC exp 1 2 t 2 C . Thus, we haved 2 (F T (Υ) t , F T (Υ ′ ) t ) ≤ C(t) t 0 d 2 (Υ s , Υ ′ s ) ds.Then, for F k T denoting k successive composition of F T , one can inductively show that it satisfiesd 2 (F k T (Υ) t , F k T (Υ ′ ) t ) ≤ (tC(t)) k k! sup s∈[0,T ] d 2 (Υ s , Υ ′ s ).Taking the supremum, we havesup s∈[0,T ] d 2 (F k T (Υ) s , F k T (Υ ′ ) s ) ≤ (T C(T )) k k! sup s∈[0,T ] d 2 (Υ s , Υ ′ s ).Thus, for a large enough k, we have shown that F k T is a contraction and from Banach Fixed Point Theorem and the completeness of the space (C([0, T ], Θ × P(R dz )), sup d), we have existence and uniqueness.Recall, the process defined in Prop. 9:dθ γ,M t = -∇ θ E γ λ (θ γ,M t , r γ,M t ) dt, where r γ,M t = 1 M M m=1 δ Z γ,M t,m ∀m ∈ [M ] : dZ γ,M t,m = b γ (θ γ,M t , r γ,M t , Z γ,M t,m) dt + 2λ r dW t,m . and γ-PVI (defined in Eq. ( 17)) augmented with extra particles (in the sense that there are M independent copies of the Z-process) to facilitate a synchronous coupling argumentdθ γ t = -∇ θ E γ λ (θ γ t , Law(Z γ t,1 )) dt, ∀m ∈ [M ] : dZ γ t,m = b γ (θ γ t , Law(Z γ t,1 ), Z γ t,m ) dt + 2λ r dW t,m .Proof of Prop. 9. This is equivalent to proving thatEQUATION)We shall treat the two terms individually. We begin with (a) in ( 24), where Jensen's inequality gives:(a) in ( 24)EQUATION)whereC θ := 2T K 2 E γ λ, we apply Cauchy-Schwarz; and the C r inequality with the Lipschitz continuity of ∇ θ E γ λ from Prop. 12. Using the C r inequality again, together with the triangle inequality:EQUATION)where rγ s = 1 M M m=1 δ Z γ s,m and we use Fournier and Guillin (2015) . Note that we also haveEQUATIONEQUATION)Applying Eq. ( 26) in Eq. ( 25) then Eqs. ( 27) and ( 28), we obtain(a) ≤ 2C θ T 0 E sup s∈[0,T ] ∥θ γ s -θ γ,M s ∥ 2 + E sup s∈[0,T ] 1 M M m=1 ∥Z γ s,m -Z γ,M s,m ∥ 2 ds + o(1). (29)Similarly, for (b) in ( 24), we have(b) = E sup t∈[0,T ] 1 M M m=1 t 0 b γ (θ γ,M s , r γ,M s , Z γ,M s,m ) -b γ (θ γ s , Law(Z γ s,1 ), Z γ s,m ) ds 2 ≤ C z E T 0 ∥θ γ,M s -θ γ s ∥ 2 + W 2 2 (r γ,M s , Law(Z γ s,1 )) + 1 M M m=1 ∥Z γ s,m -Z γ,M s,m ∥ 2 ds,where C z := 3K 2 b γ and, as before, we apply Cauchy-Schwarz, Lipschitz and C r inequalities. Then from Eqs. ( 26) to (28), we obtainEQUATION)Combining Eqs. ( 29) and ( 30) and applying Grönwall's inequality, we obtainE sup t∈[0,T ] ∥θ γ t -θ γ,M t ∥ 2 + E sup t∈[0,T ] 1 M M m=1 ∥Z γ t,m -Z γ,M t,m ∥ 2 = o(1).Taking the limit, we have the desired result.F.4 The drift in Eq. ( 17) is LipschitzIn this section, we show that the drift in the γ-PVI flow in Eq. ( 17) is Lipschitz.Proposition 11. Under the same assumptions as Prop. 8; the drift b(A, r) is Lipschitz, i.e., there exists a constant K b ∈ R >0 such that:∥ bγ (θ, r, z)-bγ (θ ′ , r ′ , z ′ )∥ ≤ K b(∥(θ, z) -(θ ′ , z ′ )∥+W 2 (r, r ′ )), ∀θ, θ ′ ∈ Θ, z, z ′ ∈ Z, r, r ′ ∈ P(Z).Proof. From definition using concavity• (which that any a, b ≥ 0, √ a + b ≤ √ a + √ b), we obtain ∥ bγ (θ, r, z) -bγ (θ ′ , r ′ , z ′ )∥ ≤ ∥∇ θ E γ λ (θ, r) -∇ θ E γ λ (θ ′ , r ′ )∥ + ∥b γ (θ, r, z) -b γ (θ ′ , r ′ , z ′ )∥. It is established below in Prop. 12 that ∇ θ E γλ satisfies a Lipschitz inequality, i.e., there is someK E γ λ ∈ R >0 such that ∥∇ θ E γ λ (θ, r) -∇ θ E γ λ (θ ′ , r ′ )∥ ≤ K E γ λ (∥θ -θ ′ ∥ + W 2 (r, r ′ )).It is established below in Prop. 13 that b γ satisfies a Lipschitz inequality, i.e., there is someK b γ ∈ R >0 such that ∥b γ (θ, r, z) -b γ (θ ′ , r ′ , z ′ )∥ ≤ K b γ (∥(θ, z) -(θ ′ , z ′ )∥ + W 2 (r, r ′ )).Hence, we have obtained as desired with K b = K E γ λ + K b γ . Proposition 12. Under the same assumptions as Prop. 8, the function (θ, r) → ∇ θ E γ λ (θ, r) is Lipschitz, i.e., there exist some constantK E γ λ ∈ R >0 such that ∥∇ θ E γ λ (θ, r) -∇ θ E γ λ (θ ′ , r ′ )∥ ≤ K E γ λ (∥θ -θ ′ ∥ + W 2 (r, r ′ )), ∀(θ, r), (θ ′ , r ′ ) ∈ M.Proof. From the definition, we have∇ θ E γ λ (θ, r) = ∇ θ E γ (θ, r) + ∇ θ R λ (θ, r). Thus, if both ∇ θ E γ and ∇ θ R λ are Lipschitz, then so is ∇ θ E γ λ .Since R λ has Lipschitz gradient (by Assumption 1), it remains to be shown that ∇ θ E γ is Lipschitz. From Prop. 6, we have∇ θ E γ (θ, r) = E p k (ϵ)r(z) (∇ θ ϕ θ • [s γ θ,r -s p ])(z, ϵ) = ∇ θ ϕ θ • d p,γ θ,r (z, ϵ) p k (dϵ)r(dz),where d p,γ θ,r (z, ϵ) := s γ θ,r (z, ϵ) -s p (z, ϵ). Then, applying Jensen's inequality, we obtainEQUATION)Focusing on the integrand, we can upper-bound it with∇ θ ϕ θ • d p,γ θ,r (z, ϵ)r(z) -∇ θ ϕ θ ′ • d p,γ θ ′ ,r ′ (z, ϵ)r ′ (z) dz (a) ≤ ∇ θ ϕ θ • d p,γ θ,r (z, ϵ) [r(z) -r ′ (z)]dz + ∇ θ ϕ θ • d p,γ θ,r (z, ϵ) -∇ θ ϕ θ ′ • d p,γ θ ′ ,r ′ (z, ϵ) r ′ (z)dz (b) ≤ ∇ θ ϕ θ • d p,γ θ,r (z, ϵ) |r(z) -r ′ (z)|dz + ∇ θ ϕ θ • d p,γ θ,r (z, ϵ) -∇ θ ϕ θ ′ • d p,γ θ ′ ,r ′ (z, ϵ) r ′ (z)dz.where in (a) we add and subtract the relevant terms and invoke the triangle inequality, and in (b) we apply Jensen's inequality. Plugging this back into Eq. ( 31), we obtainEQUATION)where the interchange of integrals is justified from Fubini's theorem for non-negative functions (also known as Tonelli's theorem).As we shall later show, the two terms have the following upper bounds:EQUATIONEQUATIONwhere K denotes a generic constant; and, upon noting that W 1 ≤ W 2 , we obtained the desired result. Now, we shall verify Eqs. ( 34) and (35). For the Eq. ( 34), we use the fact that the map z →E p k (ϵ) ∇ θ ϕ θ • d p θ,r (z, ϵ) is Lipschitz then from the dual representation of W 1 , we obtain the desired result. To see that the aforementioned map is Lipschitz, EQUATION≤ E p k (ϵ) [(a ϕ ∥ϵ∥ + b ϕ )(a d ∥ϵ∥ + b d )] ∥z -z ′ ∥ + 1 2 E p k (ϵ) (a ϕ ∥ϵ∥ + b ϕ ) 2 + d p,γ θ,r (z ′ , ϵ) 2 ∥z -z ′ ∥,EQUATIONwhich holds from the assumption that p k has finite second moments Assumption 3, and from our assumption that E p k (ϵ) d p,γ θ,r (z ′ , ϵ) is bounded. Hence, the map is Lipschitz and so Eq. ( 34) holds.As for Eq. ( 35), we focus on the integrand in Eq. ( 33)E p k (ϵ) ∇ θ ϕ θ • d p,γ θ,r (z, ϵ) -∇ θ ϕ θ ′ • d p,γ θ ′ ,r ′ (z, ϵ) ≤E p k (ϵ) ∇ θ ϕ θ • (d p,γ θ,r -d p,γ θ ′ ,r ′ )(z, ϵ) + (∇ θ ϕ θ -∇ θ ϕ θ ′ ) • d p θ ′ ,r ′ (z, ϵ) ≤E p k (ϵ) ∥∇ θ ϕ θ (z, ϵ)∥ F (d p,γ θ,r -d p,γ θ ′ ,r ′ )(z, ϵ) + ∥(∇ θ ϕ θ -∇ θ ϕ θ ′ )(z, ϵ)∥ F d p,γ θ ′ ,r ′ (z, ϵ) ≤E p k (ϵ) [(a ϕ ∥ϵ∥ + b ϕ )(a d ∥ϵ∥ + b d )] (∥θ -θ ′ ∥ + W 1 (r, r ′ )) +E p k (ϵ) (a ϕ ∥ϵ∥ + b ϕ ) d p,γ θ ′ ,r ′ (z, ϵ) ∥θ -θ ′ ∥,where, for the last line, we apply Prop. 15 and Assumption 3. Applying Young's inequality and (36), we have the desired result.Proposition 13 (b γ is Lipschitz). Under the same assumptions as Prop. 8, the map b γ is K b γ -Lipschitz, i.e., there exists a constant K b γ ∈ R >0 such that the following inequality holds for all (θ, z, r),(θ ′ , z ′ , r ′ ) ∈ Θ × R dz × P(R dz ): ∥b γ (θ, r, z) -b γ (θ ′ , r ′ , z ′ )∥ ≤ K b γ (∥(θ, z) -(θ ′ , z ′ )∥ + W 1 (r, r ′ )).Proof. One can write the drift b γ as follows (can be found in Eq. ( 43) similarly to Prop. 6), we haveb γ (θ, r, z) = -E p k (ϵ) (∇ z ϕ θ • [s γ θ,r -s p + Γ γ θ,r ])(z, ϵ) + ∇ x log p 0 (z),where Γ γ θ,r (z, ϵ) := γ∇xq θ,r (ϕ θ (z,ϵ))(q θ,r (ϕ θ (z,ϵ))+γ) 2 . Hence,∥b γ (θ, r, z) -b γ (θ ′ , r ′ , z ′ )∥ ≤ ∥E p k (ϵ) [(∇ z ϕ θ • [d p,γ θ,r + Γ γ θ,r ])(z, ϵ) -(∇ z ϕ θ ′ • [d p θ ′ ,r ′ + Γ γ θ ′ ,r ′ ])(z ′ , ϵ)]∥ + ∥∇ z log p 0 (z) -∇ z log p 0 (z ′ )∥ ≤ E p k (ϵ) ∥(∇ z ϕ θ • d p,γ θ,r + Γ γ θ,r )(z, ϵ) -(∇ z ϕ θ ′ • [d p,γ θ ′ ,r ′ + Γ γ θ ′ ,r ′ ])(z ′ , ϵ)∥ + K p0 ∥z -z ′ ∥,where for the last inequality we use Jensen's inequality and Assumption 1. Since we haveEQUATIONwhere, for (a), we add and subtract the relevant terms and invoke the triangle inequality, in (b) we use properties of the matrix norm, and in (c) we use the bounded gradient and Lipschitz gradient in Assumption 3. For Eq. (37); upon using Props. 14 and 15, which are established below, we obtainEQUATIONAs for the second term, Eq. ( 38),EQUATIONwhere (a) follows from the triangle inequality, (b) we use Prop. 14 boundedness of Γ, (c) we apply Young's inequality to the first term. Similarly to Eq. ( 36), from our Assumption 3 and our boundness assumption of the score, we have as desired. Combining Eq. ( 39) with the result of plugging Eq. ( 40) into Eq. ( 38), we obtain the result.Proposition 14 (Γ is Lipschitz and bounded). Under Assumption 2, the map (θ, r, z) → Γ γ θ,r (z, ϵ) is Lipschitz and bounded. (Lipschitz) there is constants a Γ , b Γ ∈ R >0 such that following hold:EQUATION, where x := ϕ θ (z, ϵ), and x ′ := ϕ θ ′ (z ′ , ϵ), we have:∥Γ γ θ,r (z, ϵ) -Γ γ θ ′ ,r ′ (z ′ , ϵ)∥ ≤ γ (q θ ′ ,r ′ (x ′ ) + γ)∇ x log(q θ,r (x) + γ) -(q θ,r (x) + γ)∇ x log(q θ ′ ,r ′ (x ′ ) + γ) (q θ,r (x) + γ)(q θ ′ ,r ′ (x ′ ) + γ) ≤ 1 γ |q θ ′ ,r ′ (x ′ ) -q θ,r (x)|∥∇ x log(q θ,r (x) + γ)∥ + 1 γ (q θ,r (x) + γ)∥∇ x log(q θ,r (x) + γ) -∇ x log(q θ ′ ,r ′ (x ′ ) + γ)∥ ≤ B k γ 2 |q θ ′ ,r ′ (x ′ ) -q θ,r (x)| + (B k + γ) γ ∥s γ θ,r (z, ϵ) -s γ θ ′ ,r ′ (z ′ , ϵ)∥ ≤ B k K q γ 2 (1 + a ϕ ∥ϵ∥ + b ϕ )(∥(θ, z) -(θ ′ , z ′ )∥ + W 1 (r, r ′ )) + B k + γ γ 2 (a s ∥ϵ∥ + b s )(∥(θ, z) -(θ ′ , z ′ )∥ + W 1 (r, r ′ )), where the last inequality follows from applying Prop. 18 and Assumption 3 to the first term and Prop. 16 to the last term . Hence, we have as desired.Boundedness follows from the fact thatγ∇xq θ,r (ϕ θ (z,ϵ)) (q θ,r (ϕ θ (z,ϵ))+γ) 2 ≤ 1 γ ∥∇ x q θ,r (ϕ θ (z, ϵ))∥ ≤ B k γ .Proposition 15. Under Assumptions 1 to 3, the map(θ, r, z) → s γ θ,r (z, ϵ) -s p (z, ϵ) =: d p,γ θ,r (z, ϵ) satisfies the following: there exist a d , b d ∈ R >0 such that for all (θ, r), (θ ′ , r ′ ) ∈ M, and z, z ′ ∈ R dz , we have ∥d p,γ θ,r (z, ϵ) -d p,γ θ ′ ,r ′ (z ′ , ϵ)∥ ≤ (a d ∥ϵ∥ + b d )(∥(θ, z) -(θ ′ , z ′ )∥ + W 1 (r, r ′ )).Proof. Let x := ϕ θ (z, ϵ), and x ′ := ϕ θ ′ (z ′ , ϵ). Then, we have∥d p,γ θ,r (z, ϵ) -d p,γ θ ′ ,r ′ (z ′ , ϵ)∥ ≤ ∥∇ x log p(x, y) -∇ x log p(x ′ , y)∥ + ∥∇ x log(q θ,r (x) + γ) -∇ x log(q θ ′ ,r ′ (x ′ ) + γ)∥ ≤ K p ∥x -x ′ ∥ + ∥∇ x log(q θ,r (x) + γ) -∇ x log(q θ ′ ,r ′ (x ′ ) + γ)∥ ≤ K p (a ϕ ∥ϵ∥ + b ϕ )∥(θ, z) -(θ ′ , z ′ )∥ + (a s ∥ϵ∥ + b s )(∥(θ, z) -(θ ′ , z ′ )∥ + W 1 (r, r ′ )), where Prop. 16 and Assumptions 1 and 3 are used. Proposition 16 (s is Lipschitz). Under Assumptions 2 and 3 and γ > 0, the map (θ, r, z) → s γ θ,r (z, ϵ) satisfies the following: there exist constants a s , b s ∈ R >0 such that the following inequality holds for all (θ, r), (θ ′ , r ′ ) ∈ M, and z, z ′ ∈ R dz :∥s γ θ,r (z, ϵ) -s γ θ ′ ,r ′ (z ′ , ϵ)∥ ≤ (a s ∥ϵ∥ + b s )(∥(θ, z) -(θ ′ , z ′ )∥ + W 1 (r, r ′ )),Proof. For brevity, we write x = ϕ θ (z, ϵ) and x ′ = ϕ θ ′ (z ′ , ϵ); from the definition, we haveEQUATIONwhere (a) we add and subtract the relevant terms and the triangle inequality; (b) we use the fact that∥∇ x q θ ′ ,r ′ (x ′ ) ∥ = ∥ ∇ x k θ ′ (x ′ |z)r ′ (dz)∥ ≤ B k (fromCauchy-Schwartz and the boundedness part of Assumption 2). Now, we deal with the terms individually. For the first term in Eq. ( 41), we use the fact that the map (θ, r, z) → ∇ x q θ,r (ϕ θ (z, ϵ)) is K q -Lipschitz from Prop. 17. As for the second term in Eq. ( 41), we apply Prop. 18.Hence, we obtain∥s γ θ,r (z, ϵ) -s γ θ ′ ,r ′ (z ′ , ϵ)∥ ≤ K gq γ + B k K q γ 2 (∥(θ, x) -(θ ′ , x ′ )∥ + W 1 (r, r ′ )) ≤ K gq γ + B k K q γ 2 (1 + a ϕ ∥ϵ∥ + b ϕ )(∥(θ, z) -(θ ′ , z ′ )∥ + W 1 (r, r ′ )),where we use Assumption 3 for the last line.Proposition 17. Under Assumption 2, the map (θ, r, x) → ∇ x q θ,r (x) is Lipschitz, i.e., for all ϵ, there exists a K gq ∈ R >0 such that the following inequality holds for all (θ, r), (θ ′ , r ′ ) ∈ M and z, z ′ ∈ R dz ,∥∇ x q θ,r (x) -∇ x q θ ′ ,r ′ (x ′ )∥ ≤ K gq (∥(θ, x) -(θ ′ , x ′ )∥ + W 1 (r, r ′ )).Proof. From direct computation,∥∇ x q θ,r (x) -∇ x q θ ′ ,r ′ (x ′ )∥ = [∇ x k θ (x|z)r (z) -∇ x k θ ′ (x ′ |z)r ′ (z)] dz (a) ≤ ∥∇ x [k θ (x|z) -k θ ′ (x ′ |z)]∥ r (z) dz + ∥∇ x k θ ′ (x ′ |z)∥ |r -r ′ | (z) dz (b) ≤ K gq (∥(θ, x) -(θ ′ , x ′ )∥ + W 1 (r, r ′ )),where (a) we add and subtract the appropriate terms, apply the triangle inequality and the Cauchy-Schwarz inequality; (b) for the first term, we use the Lipschitz gradient Assumption 2; and for the second term, we use the dual representation of W 1 with the fact map z → ∥∇ x log k θ (x|z)∥ is K k -Lipschitz from the reverse triangle inequality and the Lipschitz Assumption 2.Proposition 18. Under Assumption 2, the map (θ, r, x) → q θ,r (x) is Lipschitz, i.e., there exists some K q ∈ R >0 such that for all (θ, r, x), (θ ′ , r ′ , x ′ ) ∈ Θ × P(R dz ) × R dx , we have|q θ,r (x) -q θ ′ ,r ′ (x ′ )| < K q (∥(θ, x) -(θ ′ , x ′ )∥ + W 1 (r, r ′ )).Proof. From direct computation, we haveEQUATIONwhere (a) for the first term, we use the fact that the map z → |k θ (x|z)| is B k -Lipschitz (from the bounded gradient of Assumption 2), and again the Lipschitz property of k from the same assumption.G.1 Gradient estimator Proof of Prop. 6. We show the derivation of the estimators in Eq. ( 10). Eq. ( 9) will follow similarly.We have∇ z δ r E[θ, r](z) = ∇ z E k θ (x|z) log q θ,r (x) p(y, x) = ∇ z E ϵ∼p k log q θ,r (ϕ θ (z, ϵ)) p(y, ϕ θ (z, ϵ)) .Assuming that ϕ θ and p k are sufficiently regular to justify the interchange of differentiation and integration, we obtain∇ z δ r E[θ, r](z) = E ϵ∼p k ∇ z log q θ,r (ϕ(z, ϵ)) p(y, ϕ(z, ϵ)) .To obtain as desired, one can apply the chain rule.Similarly, one can derive a Monte Carlo gradient estimator for ∇ z δ r E γ as follows:∇ z δ r E γ [θ, r](z) = ∇ z E k θ (x|z) log q θ,r (x) + γ p(y, x) + q θ,r (x) q θ,r (x) + γ = ∇ z E ϵ∼p k log q θ,r (ϕ θ (z, ϵ)) + γ p(y, ϕ θ (z, ϵ)) + q θ,r (ϕ θ (z, ϵ)) q θ,r (ϕ θ (z, ϵ)) + γ .As before, if ϕ θ and p k are sufficiently regular, we obtainEQUATION)To obtain as desired, one can apply chain rule.Recall that the preconditioned gradient flow is given bydθ t = -Ψ θ t ∇ θ E λ (θ t , r t ) dt, ∂ t r t = ∇ z • (r t Ψ r t ∇ z δ r E λ [θ t , r t ]), where δ r E λ [θ, r] = δ r E[θ, r] + log r/p 0 We can rewrite the dynamics of r t as∂ t r t = ∇ z • (r t Ψ r t ∇ z [δ r E[θ t , r t ] -log p 0 + log r t ]), = ∇ z • (r t Ψ r t ∇ z [δ r E[θ t , r t ] -log p 0 ]) + ∇ z • (r t Ψ r t ∇ z log r t ).The second term can be written as ∇ z • (r t Ψ r t ∇ z log r t ) = ∇ z • (Ψ r t ∇ z r t ) = ∇ z • (∇ z • [Ψ r t r t ]) -∇ z • (r t ∇ z • Ψ r t )∂ zi    dz j=1 (Ψ r t ) ij ∂ zj r t    = dz i=1 ∂ zi    dz j=1 ∂ zj [(Ψ r t ) ij r t ] -r t ∂ zj [(Ψ r t ) ij ]    .Hence, we have the following dynamics of r t :∂ t r t = ∇ z • (r t (Ψ r t ∇ z [δ r E[θ t , r t ] -log p 0 ] -∇ z • Ψ r t )) + ∇ z • (∇ z • [Ψ r t r t ]))Examples. Following in the essence of RMSProp (Tieleman and Hinton, 2012) , we utilize the preconditioner defined as follows:B k = βB k-1 + (1 -β)Diag(A({∇ z δ r E[θ k , r k ](Z m ) 2 } M m=1 )) Ψ r k (Z) = (B k ) -0.5 where B k ∈ R dz×dz and A is some aggregation function such as the mean or max. The idea is to normalize by the aggregated gradient of the first variation across all the particles since this is the dominant component in the drift of PVI. Similarly to RMSProp, it keeps an exponential moving average of the squared gradient which can then be used in the preconditioner.In this section, we highlight additional details for reproducibility and computation. The code was written in JAX (Bradbury et al., 2018) and executed on a NVIDIA GeForce RTX 4090.H.1 Section 5.1 Hyperparameters. For the neural network, we use f θ = NN(2, 128, 2) defined in Table 3 , the number of particles M = 100, d z = 2, K = 1000, h θ = 10 -4 , h z = 10 -2 , λ r = 10 -8 , for Ψ θ we use RMSProp and we set Ψ r = I 2 .Computation Time. Each run took 8 seconds using JIT compilation.In this section, we outline all the experimental details regarding Section 5.2.Densities. Table 4 shows the densities used in the toy experiments.Banana N (z 2 ; z 2 1 /4, 1)N (z 1 ; 0, 2) X-Shape Hyperparameters. We set the number of parameter updates and particle steps to be K = 15000, and d z = 2.1 2 N 0, 2 1.8 1.8 2 + 1 2 N 0, 2 -1.8 -1.8 2 Multimodal 1 8 N 2 2 , I + 1 8 N -2 -2 , I + 1 2 N 2 -2 , I + 1 4 N -2 2 , I• f θ . We use f θ = NN(2, 512, 2).• PVI. We use M = 100, λ θ = 0, λ r = 10 -8 , h x = 10 -2 , h θ = 10 -4 , Ψ θ we use the RMSProp preconditoner, Ψ r = I dz , and L = 250.• SVI. We use K = 50 to estimate the objective (Yin and Zhou, 2018, see Eq. 5) which are around the values used in (Yin and Zhou, 2018) . We utilize RMSProp with step size 10 -4 , and r = N (0, I dz ). The implicit distribution is set to r = N (0, I dz ).• UVI. For the HMC sampler, we follow in (Titsias and Ruiz, 2019) and use 50 burn-in steps, with step-size 10 -1 and 5 leap-frog steps. We use the RMSProp optimizer with stepsize 10 -4 for k θ . The implicit distribution is set to r = N (0, I dz ).• SM. For the "dual" function written as f in the original paper (Yu and Zhang, 2023 , see Algorithm 1) we use NN(2, 512, 2). We utilize RMSProp with decaying learning rate from 10 -4 to 10 -5 to optimize the kernel k θ , and RMSProp with 10 -3 to 10 -4 for the dual function f . The implicit distribution is set to r = N (0, I dz ).Sliced Wasserstein Distance. We report the average sliced Wasserstein distance using 100 projections computed from 10000 samples from the target and the variational distribution. In this section, we outline all the hyperparameters for each method used.Hyperparameters. We use K = 20000 set d z = 10. For all kernel parameters, we use RMSProp preconditioner with step size h θ = 10 -3 .• f θ . We use f θ = NN(d z , 512, 22).• PVI. We use M = 100, λ θ = 0, λ r = 10 -8 , h x = 10 -2 , and for Ψ r we use the one described in App. G.2 with mean as the aggregate function.• SVI. We use K = 50 to estimate the objective (Yin and Zhou, 2018, see Eq. 5) which are around the values used in (Yin and Zhou, 2018) .• UVI. For the HMC sampler, we follow in (Titsias and Ruiz, 2019) and use 50 burn-in steps, with step-size 10 -1 and 5 leap-frog steps.• SM. We were unable to improve the performance of SM with our chosen kernel and instead used the implementation in https://github.com/longinYu/SIVISM?utm_source= catalyzex.com to obtain posterior samples with implementation details found in the code repository and in the paper (Yu and Zhang, 2023) .H.4 Section 5.4In this section, we outline all the experimental details regarding Section 5.4.Model. We consider the neural network BNN(d bnn in , d bnn h ) defined as fx (o) = W ⊤ 2 ReLU(W ⊤ 1 o + b 1 ) + b 2 where o ∈ R d bnn in , x = [vec(W 2 ), b 2 , vec(W 1 ), b 1 ] ⊤ , W 2 ∈ R d bnn h ×1 , b 2 ∈ R, W 1 ∈ R d bnn in ×d bnn h , b 1 ∈ R d bnn h . Given an input-output pair Y := {(O i , Y i )} B i=1 , the model can be defined as p(Y , x) = p(Y |x)p(x) where the likelihood is p(Y |x) = B i=1 N (Y i ; f x (O i ), 0.01 2 )and the prior is N (x; 0, 25I).Datasets. For all the datasets, we standardize by removing the mean and dividing by the standard deviation.• Protein. For the model, we use BNN(9, 30) Hyperparameters. We use K = 1500 set d z = 10. For all kernel parameters, we use RMSProp preconditioner with step size h θ = 10 -3 that decays to 10 -5 following a constant schedule that transitions every 100 parameters steps.• f θ , σ θ . We use f θ = NN(d z , 512, d x ) and σ θ = Softplus(NN(d z , 512, d x )) + 10 -8 and they share parameters except for the last layers. • PVI. We use M = 100, λ θ = 0, λ r = 10 -3 , h x = 10 -3 , and for Ψ r we use the one described in App. G.2 with mean as the aggregate function. • SVI. We use K = 50 to estimate the objective (Yin and Zhou, 2018, see Eq. 5) which are around the values used in (Yin and Zhou, 2018) . The implicit distribution is set to r = N (0, I dz ).• UVI. For the HMC sampler, we follow in (Titsias and Ruiz, 2019) and use 50 burn-in steps, with step-size 10 -1 and 5 leap-frog steps. The implicit distribution is set to r = N (0, I dz ). • SM. For the "dual" function written as f in the original paper (Yu and Zhang, 2023, see Algorithm 1) we use NN(d x , 512, d x ) and trained with RMSProp with stepsize 10 -2 . We tried a decaying learning schedule to 10-4 but found that this degraded the performance. We used ReLU activations instead as we found that using leaky ReLUs harmed performance. The implicit distribution is set to r = N (0, I dz ). Justification: We highlight the limitations in Section 6. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 


Decision: This paper proposes an algorithm called Particle Variational Inference (PVI) to fit semi-implicit variational distributions using gradient flows. Unlike previous SIVI approaches, PVI directly optimizes over the space of probability distributions for the mixing distribution. The paper provides solid theorical contributions, as well as numerous experiments, using multiple metrics and settings.

All four reviewers support acceptance of the paper. Reviewer JYYZ mentions the paper could potentially influence the area of generative modeling, as the proposed approximate inference tools could be used for training generative models.

However, one main weakness, shared by several reviewers, is the dimensionality of the considered problems (which is moderate but not high), which indicates that future research may be needed to target higher-dimensional distributions.

Given the novelty and good execution of the paper, I believe this is not necessarily a limitation, and recommend acceptance.

Reviews: review: paper_summary: paper_summary: The paper proposes PVI as a new method to conduct variational inference using the semi-implicit distribution. The method is to construct a gradient flow to minimize a regularized ELBO, which is practically implemented as the particle propagations. The empirical studies show the accuracy over density estimation and posterior predictions.
strength_weakness: strength_weakness: A major concern is about the simulations. First, it is unclear whether the comparisons between PVI and the baselines SVI, UVI and SM are fair. The accuracy of SIVI methods can often improve with an increase in the computation such as the number of samples in SVI and the MCMC steps in UVI. For a fair comparison, all methods need to be under the same computation budget. However, the current paper does not provide details of how the SVI, UVI, and SM are implemented. It would be better to include figures that show computation/time versus the accuracies for all the methods. 

Current simulations are relatively too simple and are only conducted on several 2D toy examples and UCI datasets. First, none of these simulations directly show the accuracy of the posterior inference for a complete Bayesian model. Second, all the simulations are in low-dimension. Does the particle method suffer from the curse of dimension problem? How does the number of particles scale with dimensionality to maintain accuracy?  It would be interesting to verify uncertainty estimation in high dimensions.

Last, the paper writing is not coherent in some places. For example
- The non-coercive is not defined in Proposition 2
- The paper mentions the non-coercive "is closely related to the problem of posterior collapse"; what is the relationship exactly?
- How to compute the precondition matrices in Eq 13 is not discussed
questions: questions: Does the regularization in Eq 4 introduce bias in the posterior inference?

Does the theoretical analysis in Sec. 4 provide guidance in designing the algorithm? 

How to check coercive and l.s.c. in Proposition 3?
limitations: limitations: N/A
paper_summary: paper_summary: The authors propose a method for SIVI called Particle Variational Inference (PVI) which employs a particle approximation of an Euclidean–Wasserstein gradient flow.  PVI directly optimizes the ELBO, and it makes no parametric assumption about the mixing distribution. Their empirical results demonstrate that PVI performs favorably against other SIVI methods across various tasks. The authors provide a theoretical analysis of the behavior of the gradient flow of a related free energy functional.
strength_weakness: strength_weakness: Experiments on a real-world application could strengthen the paper.
questions: questions: Line 329: what's the meaning of 'Here the ± denotes the average ...'?
limitations: limitations: N/A
paper_summary: paper_summary: The authors introduce a novel algorithmic approach to fitting semi-implicit variational approximations. This method is based on discretizing a suitable gradient flow, and the paper provides a comprehensive theoretical analysis to support it. The approach, called Particle Variational Inference (PVI), directly optimizes over the space of probability distributions for the mixing distribution rather than optimizing a specific parametric form for the mixing distribution, as in previous SIVI approaches.

Numerical examples demonstrate improvements compared to existing semi-implicit variational inference methods. The paper provides a solid theoretical foundation and a practical algorithm with beneficial properties over existing methods in the literature.
strength_weakness: strength_weakness: The following points, particularly the limited experimental validation, prevent me from giving the submission a higher rating (8/9):
- **Limited scope of numerical experiments:** The paper's major weakness lies in its experimental section, which represents the minimum set of experiments for an acceptable paper. Specifically: a. Lack of empirical demonstration of expressiveness: Despite claiming in Section 2 that PVI can learn potentially more expressive variational approximations than other SIVI methods, no experiments empirically demonstrate this advantage. b. Insufficient exploration of optimization stability: The paper misses an opportunity to address a key challenge in SIVI methods - the difficulty in getting objective functions/algorithms to converge (or meaningfully work). Experiments showing that PVI provides consistently more stable optimization than other SIVI methods or converges on complex models where other SIVI methods struggle would have significantly strengthened the paper. Additionally, comparing performance accuracy on Bayesian Neural Networks (BNNs) for different SIVI methods doesn't show meaningful advantages to this method.
- **Theoretical-practical gap:** The theoretical analysis is for a modified gradient flow (with γ > 0), but the practical algorithm uses γ = 0. However, this discrepancy is mitigated by the authors' justification that this approach is relatively common in the literature, and empirical results did not show obvious discrepancies between theory and practice.
- **Clarity in explaining the flexibility of PVI's variational approximation:** The novelty in PVI's approach - not optimizing a specific parametric form of the mixing distribution - was not immediately apparent to me from the explanations in Sections 2 and 3 upon first read. In particular, the last sentence in Section 2 (lines 103-105) is confusing. While it becomes more understandable after reading Section 3, the intuition behind why Q_{YuZ} might not reduce to other parameterizations when fit with PVI is not fully explained.
questions: questions: 1. How sensitive is PVI to the choice of kernel and number of particles?
2. Can the theoretical analysis be extended to cover the γ = 0 case used in practice?
3. How does PVI scale to higher-dimensional problems or larger datasets?
4. Are there specific types of problems where PVI is expected to significantly outperform existing methods?
5. How does the choice of preconditioner affect PVI's performance?
limitations: limitations: The authors have adequately addressed limitations.
paper_summary: paper_summary: This paper introduces Particle Variational Inference (PVI), a novel method for semi-implicit variational inference (SIVI) that employs empirical measures to optimize the mixing distribution without making parametric assumptions. Unlike existing SIVI methods that face challenges with intractable variational densities and rely on costly computational techniques, PVI directly optimizes the evidence lower bound (ELBO) using a particle approximation of an Euclidean–Wasserstein gradient flow.
strength_weakness: strength_weakness: I do not identify any significant weaknesses in the paper. However, there are a few points that require clarification, which I have addressed in the questions section.
questions: questions: 1.	Figure 1: The current caption makes it difficult to interpret the figure. It would be helpful to explicitly mention that lighter shades represent smaller μ\muμ values, while darker shades represent larger μ\muμ values.
2.	Figure 2: The second sentence of the caption is somewhat confusing. Referring to Table 1 for additional clarification might improve the reader's understanding.
3.	Minor Typo: In line 303, "mu=2" should be corrected to "$\mu=2$"
limitations: limitations: The authors have adequately discussed the limitations of their approach.

Comments: title: Response to rebuttal from authors
Thank you for the rebuttal
comment: Thank you for providing the additional numerical results and answering my questions. 

The paper is technically flawless and is highly novel. Semi-implicit VI has been around for years, and the construction of the objective plus algorithm solves what I consider an existing open problem in the area (unbiased gradient estimator of the ELBO directly).

Thanks for clarifying the existing numerical results. After an additional read, I am super satisfied with the benchmark versus the existing methods provided, and the wide variety of "metrics" used for assessment. While none of the target densities are in the "wow factor" territory, this is not necessary in my opinion. The paper is so strong, that I am sure that the focus of derivative works could be to "apply" or adapt the existing algorithm for high-dimensional and/or more difficult targets. 

Given that VI is the most popular approximate inference method at the moment, and that such approaches are heavily related to algorithms used to train quite a few of the popular classes of generative models at the moment, this paper has room to be of interest to many in the community, and I would some sort of spotlight or oral for this work. Thanks for the great submission.
Thank you for your response. I will keep my current score unchanged.
Thanks for the response. The score from this reviewer is unchanged.
The further analysis of the computation time and full Bayesian model, together with the discussion on the dimensionality, strengthen the paper and address my major concerns. Thanks the authors for providing code to facilitate results reproducing. I update my score accordingly.