NEO: Non Equilibrium Sampling on the Orbit of a Deterministic Transform
Year: 2022

Abstract:
Sampling from a complex distribution π and approximating its intractable normalizing constant Z are challenging problems. In this paper, a novel family of importance samplers (IS) and Markov chain Monte Carlo (MCMC) samplers is derived. Given an invertible map T, these schemes combine (with weights) elements from the forward and backward Orbits through points sampled from a proposal distribution ρ. The map T does not leave the target π invariant, hence the name NEO, standing for Non-Equilibrium Orbits. NEO-IS provides unbiased estimators of the normalizing constant and self-normalized IS estimators of expectations under π while NEO-MCMC combines multiple NEO-IS estimates of the normalizing constant and an iterated sampling-importance resampling mechanism to sample from π. For T chosen as a discrete-time integrator of a conformal Hamiltonian system, NEO-IS achieves state-of-the art performance on difficult benchmarks and NEO-MCMC is able to explore highly multimodal targets. Additionally, we provide detailed theoretical results for both methods. In particular, we show that NEO-MCMC is uniformly geometrically ergodic and establish explicit mixing time estimates under mild conditions. 35th Conference on Neural Information Processing Systems (NeurIPS 2021).

text: Sampling from a complex distribution π and approximating its intractable normalizing constant Z are challenging problems. In this paper, a novel family of importance samplers (IS) and Markov chain Monte Carlo (MCMC) samplers is derived. Given an invertible map T, these schemes combine (with weights) elements from the forward and backward Orbits through points sampled from a proposal distribution ρ. The map T does not leave the target π invariant, hence the name NEO, standing for Non-Equilibrium Orbits. NEO-IS provides unbiased estimators of the normalizing constant and self-normalized IS estimators of expectations under π while NEO-MCMC combines multiple NEO-IS estimates of the normalizing constant and an iterated sampling-importance resampling mechanism to sample from π. For T chosen as a discrete-time integrator of a conformal Hamiltonian system, NEO-IS achieves state-of-the art performance on difficult benchmarks and NEO-MCMC is able to explore highly multimodal targets. Additionally, we provide detailed theoretical results for both methods. In particular, we show that NEO-MCMC is uniformly geometrically ergodic and establish explicit mixing time estimates under mild conditions.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).

heading: 1 Introduction
text: Consider a target distribution of the form π(x) ∝ ρ(x)L(x) where ρ is a probability density function (pdf) on Rd and L is a nonnegative function. Typically, in a Bayesian setting, π is a posterior distribution associated with a prior distribution ρ and a likelihood function L. Another situation of interest is generative modeling where π is the distribution implicitly defined by a Generative Adversarial Networks (GAN) discriminator-generator pair where ρ is the distribution of the generator and L is derived from the discriminator [31, 6]. In a Variational Auto Encoder (VAE) context [14, 5], π could be the true posterior distribution, ρ the approximate posterior distribution output by the encoder and L an importance weight between the true posterior and approximate posterior distributions. We are interested in this paper in sampling from π and approximating its intractable normalizing constant Z = ∫ ρ(x)L(x)dx. These problems arise in many applications in statistics, molecular dynamics or machine learning, and remain challenging.
Many approaches to compute normalizing constants are based on Importance Sampling (IS) - see [1, 2] and the references therein - and its variations, among others, Annealed Importance Sampling (AIS) [19, 34, 10] and Sequential Monte Carlo (SMC) [9]. More recently, Neural IS has also become very popular in machine learning; see e.g. [11, 18, 21, 22, 32, 33]. Neural IS is an adaptive IS which relies on an importance function obtained by applying a normalizing flow to a reference distribution. The parameters of this normalizing flow are chosen by minimizing a divergence between the proposal and the target (such as the Kullback–Leibler [18] or the χ2-divergence [1]). Recent work on the subject proposes to add stochastic moves in order to enhance the performance of the normalizing flows [33].
More recently, the Non-Equilibrium IS (NEIS) method has been introduced by [23] as an alternative to these approaches. Similar to Neural IS, NEIS consists in transporting samples {Xi}Ni=1 from a reference distribution using a family of deterministic mappings. For NEIS, this family is chosen to be an homogeneous differential flow (φt)t∈R. In contrast to Neural IS, for any i ∈ [N ], the sample Xi is propagated both forward and backward in time along the orbits associated with (φt)t∈R until stopping conditions are met. Moreover, the resulting estimator of the normalizing constant is obtained by computing weighted averages of the whole orbit (φt(Xi))t∈[τ+,i,τ−,i], where τ+,i, τ−,i are the resulting stopping times, and not only the endpoints φτ+,i(X i), φτ−,i(X i). In [23], the authors provide an application of NEIS with (φt)t∈R associated with a conformal Hamiltonian dynamics, and reports impressive numerical results on difficult normalizing constants estimation problems, in particular for high-dimensional multimodal distributions.
We propose in this work NEO-IS which alleviates the shortcomings of NEIS. Similar to NEIS, samples are drawn from a reference distribution, typically set to ρ, and are propagated under the forward and backward orbits of a discrete-time dynamical system associated with an invertible transform T. An estimator of the normalizing constant is obtained by reweighting all the points on the whole orbits using the IS rule. Contrary to NEIS, the NEO-IS estimator of Z is unbiased under assumptions that are mild and easy to verify. It is more flexible than NEIS because it does not rely on the accuracy of the discretization of a continuous-time dynamical system.
We then show how it is possible to leverage the unbiased estimator of Z defined by NEO-IS to obtain NEO-MCMC, a novel massively parallel MCMC algorithm to sample from π. In a nutshell, NEOMCMC relies on parallel walkers which each estimates the normalizing constant but are allowed to interact through a resampling mechanism.
Our contributions can be summarized as follows.
(i) We present a novel class of IS estimators of the normalizing constant Z referred to as NEO-IS. More broadly, a small modification of this algorithm also allows us to estimate integrals with respect to π. Both finite sample and asymptotic guarantees are provided for these two methodologies.
(ii) We develop a new massively parallel MCMC method, NEO-MCMC. NEO-MCMC combines NEO-IS unbiased estimator of the normalizing constant with iterated samplingimportance resampling methods. We prove that it is π-reversible and ergodic under very general conditions. We derive also conditions which imply that NEO-MCMC is uniformly geometrically ergodic (with an explicit expression of the mixing time).
(iii) We illustrate our findings using numerical benchmarks which show that both NEO-IS and NEO-MCMC outperform state-of-the-art (SOTA) methods in difficult settings.

heading: 2 NEO-IS algorithm
text: In this section, we derive the NEO-IS algorithm. The two key ingredients for this algorithm are (1) the reference distribution ρ and (2) a transformation T assumed to be a C1-diffeomorphism with inverse T−1. Write, for k ∈ N∗ = N\{0}, Tk = T ◦Tk−1, T0 = Idd and similarly T−k = T−1 ◦T−(k−1). For any k ∈ Z, denote by ρk : Rd → R+ the pushforward of ρ by Tk, defined for x ∈ Rd by ρk(x) = ρ(T
−k(x))JT−k(x), where JΦ(x) ∈ R+ is the absolute value of the Jacobian determinant of Φ : Rd → Rd evaluated at x. In line with multiple importance sampling à la Owen and Zhou [20], we introduce the proposal density
ρT(x) = Ω −1 ∑
k∈Z $kρk(x) , (1) where {$k}k∈Z is a nonnegative sequence and Ω = ∑ k∈Z$k. Note that we assume in the sequel that the support of the weight sequence defined as {k ∈ Z : $k 6= 0} is finite. Thus, the mixture distribution in (1) is a finite mixture. Given x ∈ Rd, ρT(x) is a function of the forward and backward orbit of T through x.
For any nonnegative function f , the definition of ρT implies that ∫ f(y)ρT(y)dy =
Ω−1 ∫ ∑
k∈Z$kf(T k(x))ρ(x)dx. Assuming that $0 > 0, the ratio ρ(x)/ρT(x) ≤ $−10 Ω <∞ is
bounded. We can therefore apply the IS principle which allows to write the identity∫ f(x)ρ(x)dx = ∫ ( f(y) ρ(y)
ρT(y)
) ρT(y)dy = ∫ ∑ k∈Z f(Tk(x))wk(x)ρ(x)dx , (2)
where the weights are given by (see Appendix A.2 for a detailed derivation),
wk(x) = $kρ(T k(x))/{ΩρT(Tk(x))} = $kρ−k(x) /∑ i∈Z $k+iρi(x) . (3)
We assume in the sequel that $0 > 0. In particular, note that under this condition, the weights wk are also upper bounded uniformly in x: for any x ∈ Rd, wk(x) ≤ $k/$0. Equations (2) and (3) suggest to estimate the integral ∫ f(x)ρ(x)dx by INEO$,N (f) = N −1∑N i=1 ∑ k∈Z wk(X i)f(Tk(Xi)) where {Xi}Ni=1 are i.i.d. samples from the proposal ρ, which is denoted by X1:N iid∼ ρ.
Algorithm 1 NEO-IS Sampler
1. Sample X1:N iid∼ ρ for i ∈ [N ]. 2. For i ∈ [N ], compute the path (Tj(Xi))j∈Z and weights (wj(Xi))j∈Z.
3. INEO$,N (f) = N −1∑N i=1 ∑ k∈Z wk(X i)f(Tk(Xi)).
This estimator is obtained by a weighted combination of the elements of the independent forward and backward orbits {Tk(Xi)}k∈Z with X1:N
iid∼ ρ. This estimator is referred to as NEO-IS. Choosing f ≡ L provides the NEO-IS estimator of the normalizing constant of π:
Ẑ $ Xi = ∑ k∈Z L(T k(Xi))wk(X i) , Ẑ $ X1:N = N −1∑N i=1 Ẑ $ Xi . (4)
We now study the performance of the NEO-IS estimator. The following two quantities play a fundamental role in the analysis:
E$T = EX∼ρ [(∑ k∈Z wk(X)L(T k(X))/Z )2] ,M$T = supx∈Rd ∑ k∈Z wk(x)L(T k(x))/Z . (5)
Theorem 1. Ẑ $
X1:N is an unbiased estimator of Z. If E $ T < ∞, then, E[|Ẑ
$
X1:N /Z−1|2] = N−1(E$T − 1). If M$T <∞, then, for any δ ∈ (0, 1), with probability 1− δ, √ N ∣∣∣Ẑ$X1:N /Z−1∣∣∣ ≤
M$T √ log(2/δ)/2.
The (elementary) proof is postponed to Appendix A.3. E$T plays the role of the second-order moment of the importance weights EX∼ρ[L2(X)] which is key to the performance of IS algorithms
[1, 2]. In addition, since the NEO-IS estimator Ẑ $
X1:N is unbiased, the Cauchy–Schwarz inequality implies that EX∼ρ [(∑ k∈Z wk(X)L(T k(X)))2 ] ≥ Z2 and hence that E$T ≥ 1. Note that if ‖L‖∞ = supx∈Rd L(x) < ∞, then since the weights are uniformly bounded by Ω$−10 , we have M$T ≤ ‖L‖∞Ω$ −1 0 /Z.
Using the NEO-IS estimate Ẑ $
X1:N of the normalizing constant, we can construct a self-normalized IS estimate of ∫ f(x)π(x)dx:
JNEO$,N (f) = N −1 N∑ i=1 Ẑ $ Xi Ẑ $ X1:N ∑ k∈Z L(Tk(Xi))wk(X i) Ẑ $ Xi f(Tk(Xi)) , (6)
referred to as NEO-SNIS estimator. This expression may seem unnecessarily complicated but highlights the hierarchical structure of the estimator. We combine estimators (Ẑ $ Xi) −1∑ k∈Z L(T k(Xi))wk(X
i)f(Tk(Xi)) evaluated on the forward and backward orbits through the points {Xi}Ni=1 using the normalized weights {Ẑ $ Xi/Ẑ $
X1:N }Ni=1. Although the NEO-IS estimator is unbiased, the NEO-SNIS is in general biased. However, for bounded functions, both the bias and the variance of the NEO-SNIS estimator are O(N−1), with constants proportional to E$T . For g a π-integrable function, we set π(g) = ∫ g(x)π(x)dx.
Theorem 2. Assume that E$T <∞. Then, for any function g satisfying supx∈Rd |g(x)| ≤ 1 on Rd, and N ∈ N
E X1:N iid∼ρ
[ |JNEO$,N (g)− π(g)|2 ] ≤ 4 ·N−1E$T , (7)∣∣∣E
X1:N iid∼ρ
[ JNEO$,N (g)− π(g) ]∣∣∣ ≤ 2 ·N−1E$T . (8) If M$T <∞, then for δ ∈ (0, 1], with probability at least 1− δ,√
N |JNEO$,N (g)− π(g)| ≤ ‖g‖∞M$T √ 32 log(4/δ) . (9)
The proof is postponed to Appendix A.4. These results extend to NEO-SNIS estimators the results known for self-normalized IS estimators; see e.g., [1, 2] and the references therein. The upper bounds stated in this result suggest it is good practice to keep E$T /N small in order to obtain sensible approximations. For two pdfs p and q on Rd, denote by Dχ2(p, q) = ∫ {p(x)/q(x)− 1}2q(x)dx the χ2-divergence between p and q. Lemma 3. For any nonnegative sequence ($k)k∈Z, we have E$T ≤ Dχ2(π‖ρT) + 1.
The proof is postponed to Appendix A.5. Lemma 3 suggests that accurate sampling requires N to scale linearly with the χ2-divergence between the target π and the extended proposal ρT.
Remark 1. We can extend NEO to non homogeneous flows, replacing the family {Tk : k ∈ Z} with a collection of mappings {Tk : k ∈ Z}. This would allow us to consider further flexible classes of transformations such as normalizing flows; see e.g. [21]. The χ2-divergence Dχ2(π‖ρT) provides a natural criterion for learning the transformation. We leave this extension to future work.
Conformal Hamiltonian transform The efficiency of NEO relies heavily on the choice of T. Intuitively, a sensible choice of T requires that (i) E$T is small, i.e. ρT should be close to π by Lemma 3 (see (5)), (ii) the inverse T−1 and the Jacobian of T are easy to compute. Following [23], we use for T a discretization of a conformal Hamiltonian dynamics. Assume that U(·) = − log π(·) is continuously differentiable. We consider the augmented distribution π̃(q, p) ∝ exp{−U(q)−K(p)} on R2d, where q is the position, p is the momentum, and K(p) = pTM−1p/2 is the kinetic energy, with M a positive definite mass matrix. By construction, the marginal distribution of the momentum under π̃ is the target pdf π(q) = ∫ π̃(q, p)dp. The conformal Hamiltonian ODE associated with π̃ is defined by
dqt/dt = ∇pH(qt, pt) = M−1pt , (10) dpt/dt = −∇qH(qt, pt)− γpt = −∇U(qt)− γpt ,
where H(q, p) = U(q) +K(p), and γ > 0 is a damping constant. Any solution (qt, pt)t≥0 of (10) satisfies setting Ht = H(qt, pt), dHt/dt = −γpTt M−1pt ≤ 0. Hence, all orbits converge to fixed points that satisfy∇U(q) = 0 and p = 0; see e.g. [12, 17].
In the applications below, we consider the conformal version of the symplectic Euler (SE) method of (10), see [12]. This integrator can be constructed as a splitting of the two conformal and conservative parts of the system (10). When composing a dissipative with a symplectic operator, we set for all (q, p) ∈ R2d, Th(q, p) = (q + hM−1{e−hγp − h∇U(q)}, e−hγp − h∇U(q)), where h > 0 is a discretization stepsize. This transformation can be connected with classical momentum optimization schemes, see [12, Section 4]. By [12, Section 3], for any h > 0 Th is a C1-diffeomorphism on R2d with Jacobian given by JTh(q, p) = e−γhd. In addition, its inverse is T −1 h (q, p) = (q − hM−1p, eγh{p+ h∇U(q − hM−1p)}). Therefore, the weight (3) of the NEO estimator is given by
wk(q, p) = $kρ̃(T
k h(q, p))e −γkhd∑ j∈Z$k+j ρ̃(T −j h (q, p))e γjhd ,
where ρ̃(q, p) ∝ ρ(q)e−K(p). Figure 1 displays for different values of γ on a log-scale the bound E 1[0:K]
Th − 1 appearing in Theorem 1 as a function of K, here we use the sequence of weights ($k)k∈Z = (1[0:K](k))k∈Z (i.e. only the K + 1 first elements of the forward orbits are used and are equally weighted). For comparison, we also present on the same plot the bounds achieved by averaging K + 1 independent IS estimates, EIS(K)− 1 = (K + 1)−1EX∼ρ[L(X)2]. Interestingly, Figure 1 shows that there is a trade-off in the choice of γ which controls the exploration of the state space by the Hamiltonian dynamics since the higher γ, the faster the orbits converge towards the modes. This fast convergence prevents a “good” exploration of the space; e.g. E 1[0:K]
Th is smaller for
γ = 1.0 than for γ = 2.0 when K > 7.
3 NEO-MCMC algorithm
We now derive an MCMC method to sample from π based on the NEO-IS estimator. A natural idea consists in adapting the Sampling Importance Resampling procedure (SIR) (see for example [24, 27]) to the NEO framework.
The SIR method to sample JNEO$,N (see (6)) consists of 4 steps. (SIR-1) Draw independently X1:N iid∼ ρ and compute the associated forward and backward orbits {Tk(Xi)}k∈Z of the point. (SIR-2) Compute the normalizing constants associated with each orbit {Ẑ $
Xi}Ni=1. (SIR-3) Sample an orbit index IN ∈ [N ] with probability {Ẑ $ Xi/ ∑N j=1 Ẑ $
Xj}Ni=1. (SIR-4) Draw the iteration index KN on the IN -th orbit with probability {L(Tk(XIN ))wk(XI N )/Ẑ $ XIN }k∈Z.
The resulting draw is denoted by UN = TK N (XI N
). By construction, for any bounded function f , we get that E [ f(UN ) ∣∣X1:N , IN ] = {Ẑ$XIN }−1∑k∈Z wk(XIN )L(Tk(XIN )) which implies E [ f(UN )
∣∣X1:N ] = JNEO$,N (f) (see (6)). Using Theorem 2, we therefore obtain |E[f(UN )] −∫ f(z)π(z)dz| ≤ 101/2‖f‖∞E$T N−1, showing that the law of the random variable µN = Law(UN ) converges in total variation to π as N →∞,
‖µN − π‖TV = sup ‖f‖∞≤1 |µN (f)− π(f)| ≤ 101/2E$T N−1 . (11)
Algorithm 2 NEO-MCMC Sampler At step n ∈ N∗, given the conditioning orbit point Yn−1. Step 1: Update the conditioning point
1. Set X1n = Yn−1 and for any i ∈ {2, . . . , N}, sample Xin
iid∼ ρ. 2. Sample the orbit index In with probability
proportional to (Ẑ $
Xin )i∈[N ], (4).
3. Set Yn = XInn . Step 2: Output a sample
4. Sample index Kn with probability proportional to {wk(Yn)L(Tk(Yn))/Ẑ $
Yn}k∈Z 5. Output Un = TKn(Yn).
Based on [3], we now derive the NEOMCMC procedure, which in a nutshell consists in iterating the SIR procedure while keeping a conditioning point (or equivalently, orbit); see Appendix C. The convergence of NEO-MCMC does not rely on letting N →∞: the NEO-MCMC works as soon as N ≥ 2, although as we will see below the mixing time decreases as N increases.
This procedure is summarized in Algorithm 2. The NEO-MCMC procedure is an iterated algorithm which produces a sequence {(Yn, Un)}n∈N of points in Rd. The n-th iteration of the NEO-MCMC algorithm consists in two main steps: 1) updating the conditioning point Yn−1 → Yn 2) sampling Un by selecting a point in the orbit {Tk(Yn)}k∈Z of the conditioning point. Compared to SIR, only the generation of the points (step (SIR-1)) is modified: we set X1n = Yn−1 (the conditioning point), and then draw X2:Nn
iid∼ ρ. The sequence {Yn}n∈N defined by Algorithm 2 is a Markov chain: P (Yn ∈ A |Y0:n−1) = P (Yn ∈ A |Yn−1) = P (Yn,A) where
P (y,A) = ∫ δy(dx 1) N∏ j=2 ρ(xj)dxj N∑ i=1
Ẑ $ xi∑N j=1 Ẑ $ xj 1A(x i) , y ∈ Rd ,A ∈ B(Rd) . (12)
Note that this Markov kernel describes the way, at stage n + 1, the conditioning point Yn+1 is selected given Yn, which depends only on the estimator of the normalizing constants associated with each orbit, but not on the sample Un selected on the conditioning orbit. In addition, given the conditioning point Yn at the n-th iteration, the conditional distribution of the output sample Un is P ( Un ∈ B | In, X1:Nn ) = P (Un ∈ B |Yn) = Q(Yn,B) where
Q(y,B) = ∑ k∈Z wk(y)L(T k(y)) Ẑ $ y 1B(T k(y)) , y ∈ Rd ,B ∈ B(Rd) . (13)
With these notations, if the Markov chain is started at Y0 = y, then for any n ∈ N, the law of the n-th conditioning point is P (Yn ∈ A |Y0 = y) = Pn(y,A) and the law of the n-th sample is P (Un ∈ B |Y0) = PnQ(y,B). Define π̃ the pdf given for y ∈ Rd by
π̃(y) = ρ(y)
Z ∑ k∈Z wk(y)L(T k(y)) = ρ(y)Ẑ $ y Z . (14)
The following theorem shows that, for any initial condition y ∈ Rd, the distribution of the variable Yn converges in total variation to π̃ and that the distribution of Un converges to π. Theorem 4. The Markov kernel P is reversible with respect to the distribution π̃, ergodic and Harris positive, i.e., for all y ∈ Rd, limn→∞ ‖Pn(y, ·) − π̃‖TV = 0. In addition, π = π̃Q and limn→∞ ‖PnQ(y, ·) − π‖TV = 0. Moreover, for any bounded function g and any y ∈ Rd, limn→∞ n −1∑n−1 i=0 g(Ui) = π(g), P-almost surely, where {Ui}i∈N is defined in Algorithm 2 with Y0 = y.
The proof is postponed to Appendix A.6. Remark 2. We may provide another sampling procedure of {Yn}n∈N. Define the pdf on the extended space [N ]× RdN by π̌(i, x1:N ) = N−1π̃(xi) ∏N j=1,j 6=i ρ(x j). Consider a Gibbs sampler
targeting π̌ consisting in (a) sampling X1:N\{In−1}n |(In−1, Xn−1) ∼ ∏ j 6=In−1 ρ(x j), (b) sampling In|X1:Nn ∼ Cat({Ẑ $ Xin / ∑N j=1 Ẑ $ Xjn }Ni=1 and (c) set Yn = XInn . This algorithm is a Gibbs sampler on π̌ and we easily verify that the distribution of {Yn}n∈N is the same as Algorithm 2.
The next theorem provides non asymptotic quantitative bounds on the convergence in total variation. The main interest of NEO-MCMC algorithm is motivated empirically from observed behaviour: the mixing time of the corresponding Markov chain improves as N increases. This behaviour is quantified theoretically in the next theorem. Moreover, this improvement is obtained with little extra computational overhead, since sampling N points from the proposal distribution ρ, computing the forward and backward orbits of the points and evaluating the normalizing constants {Ẑ $
Xin }Ni=1 can
be performed in parallel. Theorem 5. Assume thatM$T <∞, see (5). Set N = (N − 1)/(2M$T +N − 2) and κN = 1− N . Then, for any y ∈ Rd and k ∈ N, ‖P k(y, ·)− π̃‖TV ≤ κkN and ‖P kQ(y, ·)− π‖TV ≤ κkN .
Instead of sampling the new points X2:Nn independently from ρ (Step 1 in Algorithm 2), it is possible to draw the proposals X1:Nn conditional to the current point Yn−1; see [29, 8, 26, 25] for related works. Following [25], we use a reversible Markov kernel with respect to the proposal ρ, i.e., such that ρ(x)m(x, x′) = ρ(x′)m(x′, x), assuming for simplicity that this kernel has density m(x, x′). If ρ = N(0, σ2 Idd) , an appropriate choice is an autoregressive kernel m(x, x′) = N(x′;αx, σ2(1− α2) Idd). More generally, we can use a Metropolis–Hastings kernel with invariant distribution ρ. In this case, r1(x1, x1:N\{1}) = ∏N j=2m(x j−1, xj) and for each i ∈ [2 : N ],
ri(x i, x1:N\{i}) = i−1∏ j=1 m(xj+1, xj) N∏ j=i+1 m(xj−1, xj) . (15)
Since m is reversible with respect to ρ, for all i, j ∈ [N ], ρ(xi)ri(xi, x1:N\{i}) = ρ(xj)rj(x j , x1:N\{j}) where ri(xi;x1:N\{i}) defines the the conditional distribution of X1:N\{i} given Xi = xi. The only modification in Algorithm 2 is Step 1, which is replaced by: Draw Un ∈ [N ] uniformly, set XUnn = Yn−1 and sample X 1:N\{Un} n ∼ rUn(XUnn , ·). The validity of this procedure is established in Appendix A.6.
4 Continuous-time version of NEO and NEIS
The NEO framework can be thought of as an extension of NEIS introduced in [23]. NEIS focuses on normalizing constant estimation and should be therefore compared with NEO-IS. In [23], the authors do not consider possible extensions of these ideas to sampling problems. We consider here how NEO could be adapted to continuous-time dynamical system. Proofs of the statements and detailed technical conditions are postponed to Appendix B.
Consider the Ordinary Differential Equation (ODE) ẋt = b(xt) , where b : Rd → Rd is a smooth vector field. Denote by (φt)t∈R the flow of this ODE (assumed to be well-behaved). Under appropriate regularity condition Jφt(x) = exp( ∫ t 0 ∇ · b(φs(x))ds); see Lemma S5. Let $ : R → R+ be a
nonnegative smooth function with finite support, with Ωc = ∫∞ −∞$(t)dt. The continuous-time counterpart of the proposal distribution (1) is ρcT(x) = (Ω c)−1 ∫∞ −∞$(t)ρ(φ−t(x))Jφ−t(x)dt, which is a continuous mixture of the pushforward of the proposal ρ by the flow of (φs)s∈R. Assuming for simplicity that ρ(x) > 0 for all x ∈ Rd, then ρcT(x) > 0 for all x ∈ Rd, and using again the IS formula, for any nonnegative function f ,∫
f(x)ρ(x)dx = ∫ f(x) ρ(x)
ρcT(x) ρcT(x)dx = ∫ [∫ ∞ −∞ wct (x)f(φt(x))dt ] ρ(x)dx, (16)
wct (x) = $(t)ρ(φt(x))Jφt(x) /∫ ∞ −∞ $(s+ t)ρ(φs(x))Jφs(x)ds . (17)
These relations are the continuous-time counterparts of (2). Eqs. (16)-(17) define a version of NEIS [23], with a finite support weight function $; see Appendices B.2 and B.3 for weight functions with infinite support. This identity is of theoretical interest but must be discretized to obtain a computationally tractable estimator. For h > 0, denote by Th an integrator with stepsize h > 0 of the ODE ẋ = b(x). We may construct NEO-IS and NEO-SNIS estimators based on the transform T ← Th and weights $k ← $(kh). We might show that for any bounded function f and for any x ∈ Rd, limh↓0 ∑ k∈Z wk(x)f(T k h(x)) = ∫∞ −∞ w c t (x)f(φt(x))dt, where we omitted here the
dependency in h of wk. Therefore, taking h ↓ 0+, the NEO-IS converges to the continuous-version (16)-(17). There is however an important difference between NEO and the NEIS method in [23] which stems from the way (16)-(17) are discretized. Compared to NEIS, NEO-IS using T← Th and weights $k ← $(kh) is unbiased for any stepsize h > 0. NEIS uses an approach inspired by the nested-sampling approach, which amounts to discretizing the integral in (16) also in the state-variable x; see [28, 7]. This discretization is biased which prevents the use of this approach to develop MCMC sampling algorithm; see Appendix B.

heading: 5 Experiments and Applications
text: Normalizing constant estimation The performance of NEO-IS is assessed on different normalizing constant estimation benchmarks; see [13]. We focus on two challenging examples. Additional experiments and discussion on hyperparameter choice are given in the supplementary material, see Appendix D.1.
(1) Mixture of Gaussian (MG25): π(x) = P−1 ∑P i=1 N(x;µi,j , Dd), where d ∈ {10, 20, 45}, Dd = diag(0.01, 0.01, 0.1, . . . , 0.1) and µi,j = [i, j, 0, . . . , 0]T with i, j ∈ {−2, . . . , 2}.
(2) Funnel distribution (Fun) π(x) = N(x1; 0, a2) ∏d i=1 N(xi; 0, e
2bx1) with d ∈ {10, 20, 45}, a = 1, and b = 0.5. In both case, the proposal is ρ = N(0, σ2ρ Idd) with σ 2 ρ = 5.
The NEO-IS estimator is compared with (i) the IS estimator using the proposal ρ, (ii) the Adaptive Importance Sampling (AIS) estimator of [30], (iii) Stochastic Normalizing Flows (SNF)1 and (iv) the Neural Importance Sampling (NIS)2. For NEO-IS, we use $k = 1[K](k) with K = 10 (ten steps on the forward orbit), and conformal Hamiltonian dynamics γ = 1, M = 5 · Idd for dimensions d = {10, 20}, and γ = 2.5 for d = 45 (where γ is the damping factor, M the mass matrix, h is the stepsize of the integrator). The parameters of AIS are set to obtain a complexity comparable to NEO-IS; see Appendix D.1. For NIS, we use the default parameters and for SNF we used the same architectures as in [33]. In Fun, we set γ = 0.2, K = 10, M = 5 · Idd, and h = 0.3. The IS estimator was based on 5 · 105 samples, and NIS, NEO-IS and AIS were computed with 5 · 104 samples. Figure 2 shows that NEO-IS consistently outperforms the competing methods.
1Implementation available at https://github.com/noegroup/stochastic_normalizing_flows. 2Implementation available at https://github.com/ndeutschmann/zunis.
Sampling NEO-MCMC is assessed for the distributions (MG25) (d = 40) and Fun (d = 20). NEO-MCMC sampler is compared with (i) the No-U-Turn Sampler - Pyro library [4] - and (ii) i-SIR algorithm [25]. The proposal distribution is ρ = N(0, σ2ρ Idd) with σ 2 ρ = 5. Dependent proposals are used (see (15)) with m(x, x′) = N(x′;αx, σ2ρ(1− α2) Idd) with α = 0.99. For NUTS, the default parameters are used. For i-SIR, we use the same number of proposals N = 10, proposal distribution and dependent proposal as for NEO-MCMC. To perform a fair comparison, we use the same clock time for all three algorithms. The number of iterations for correlated i-SIR, NEO-MCMC, and NUTS are n = 4 · 106, n = 4 · 105, and n = 5 · 105, respectively. Figure 3 displays the empirical
two-dimensional histograms of the two first coordinates of samples from the ground truth, i-SIR, NUTS and NEO-MCMC sampler. It is worthwhile to note that NEO-MCMC algorithm performs much better for MG25 which is a very challenging distribution, even for SOTA algorithm such as NUTS, which struggles to cross energy barriers between modes. For Fun, NEO-MCMC performs favourably with respect to NUTS, which is well adapted for this type of distributions.
Block Gibbs Inpainting with Deep Generative models and NEO-MCMC We apply NEOMCMC to the task of sampling the posterior of a deep latent variable model. To be consistent with the rest of the paper, we use non-standard notation here with x being the latent variable and z the observation. More precisely, we assume that x ∼ N(0, Idd) and a conditional distribution p(z | x) which generates an image z = (z1, . . . , zD) ∈ RD. Given a family of parametric decoders {x 7→ pθ(z | x), θ ∈ Θ}, and a training set D = {zi}Mi=1, training involves finding the MLE θ∗ = arg maxθ∈Θ pθ(D). As pθ(z) = ∫ pθ(z | x)p(x)dx, the likelihood is intractable and to alleviate this problem, [14] proposed to train jointly an approximate posterior qφ(x|z) that maximizes a tractable lower-bound on the log-likelihood: ELBO(z, θ, φ) = EX∼qφ(·|z)[log pθ(z,X)/qφ(X|z)] ≤ pθ(z), where qφ(x | z) is a tractable conditional distribution with parameters φ ∈ Φ. It is assumed in the sequel that conditional to the latent variable x, the coordinates are independent, i.e. pθ(z | x) = ∏D i=1 pθ(z
i|x). Note that it is possible to train VAE with the NEO algorithm, using the unbiased estimate of the normalizing constant to construct an ELBO. This approach is described in the supplement Appendix E. We do not focus on this approach here and assume that the VAE has been trained and we are only interested in the sampling problem. In our experiment, we use a VAE trained on CelebA dataset 3 [16]. We consider the Block Gibbs inpainting task introduced in [15, Section 5.2.2]. Given an image z, denote by [zt, zb] the top and the bottom half pixels. Assume only zt? is observed, then we are interested in in-painting the bottom of an image by the posterior distribution of zb given zt?. This is achieved using Block Gibbs sampling. A two-stage Gibbs sampler amounts to (a) sampling pθ∗(x|zt, zb) and (b) sampling pθ∗(zb|x, zt) = pθ∗(zb|x) (since zb and zt are independent conditional on x). Given zk = (zt?, z b k), we sample at each step xk ∼ pθ∗(x | zk) and then
3Publicly available online, see https://github.com/YannDubs/disentangling-vae/tree/master/ results/betaH_celeba
zbk+1 ∼ pθ∗(zb | xk). We then set zk+1 = (zt∗, zbk+1). Stage (b) is elementary but stage (a) is challenging. We use an MCMC-within-Gibbs scheme using different samplers. We use the following decomposition of pθ∗(x | z) ∝ ρ(x)L(x) for ρ(x) ∝ qβφ∗(x | z) and L(x) = pθ∗(x, z)/q β φ∗(x | z) with β ∈ (0, 1). It is possible to sample from ρ(x) as qφ∗(x | z) is Gaussian. In our experiments with CelebA and the chosen trained VAE, we have x ∈ R10 (recall that x is our latent variable here), z ∈ R12288, and use β = 0.1. We then compare i-SIR, HMC and NEO-MCMC sampler in stage (a), with the same computational complexity (N = 10, K = 12, γ = 0.2 for NEO-MCMC, N = 120 for i-SIR, and HMC is run with K = 20 leap-frog steps). Again, NEO-MCMC and i-SIR use dependent proposals, with m a Random Walk Metropolis kernel with stepsize 0.1. For each algorithm, 10 steps are performed. Figure 4 displays the evolution of the resulting Markov chains. The samples clearly illustrate that NEO-MCMC mixes better than i-SIR and HMC. More details and examples are presented in the supplementary.

heading: 6 Conclusion
text: In this paper, we have proposed a new family of algorithms, NEO, for computing normalizing constants and sampling from complex distributions. This methodology comes with asymptotic and non-asymptotic convergence guarantees. For normalizing constant estimation, NEO-IS compares favorably to state-of-the-art algorithms on difficult benchmarks. NEO-MCMC is also able to sample some complex distributions: it is particularly well-adapted to sampling multimodal distributions, thanks to its proposal mechanism which avoids being trapped in local modes. There are numerous potential extensions to this work. For example, it would be interesting to consider deterministic transformations other than conformal Hamiltonian dynamics integrators. These transformations could be trained, as for Neural IS, using a variation lower bound. It would also be interesting to further investigate the influence of the mixture weights {$k}k∈Z on the efficiency of NEO. Broader impact: Sampling from complex target distributions and computing their normalizing constants has numerous applications. Our work proposes novel methods to address such problems and has thus potential applications in many areas. This work does not present any foreseeable societal consequence.
Acknowledgments and Disclosure of Funding
Arnaud Doucet is partly supported by the EPSRC grant EP/R034710/1. He also acknowledges support of the UK Defence Science and Technology Laboratory (DSTL) and EPSRC under grant EP/R013616/1. This is part of the collaboration between US DOD, UK MOD and UK EPSRC under the Multidisciplinary University Research Initiative. Alain Durmus and Eric Moulines acknowledge support of the Lagrange Mathematical and Computing Research Center.

Review summaries:
Based on Non-Equilibrium IS (NEIS), the authors propose the new importance sampler NEO-IS. The main advantage of NEO-IS over NEIS seems to be that NEO-IS provides unbiased estimates of the normalization constant. They prove the unbiasedness and also propose an MCMC extension for sampling (which was not provided by NEIS). They empirically show that NEO-IS (and its MCMC variant) can be better that other recent methods for tasks of  (1) finding the marginal likelihood, (2) Sampling from posterior distribution, (3) Inpaiting. 

The paper proposes a new family of importance samplers based on iterating an invertible map $T$ forwards and backwards, thus generating orbits that serve as importance samples after proper weighting. If normalized importance weights can be computed, the resulting IS estimate is unbiased. Also a self-normalized version is proposed, which is biased but still converges against the correct expectation for $N\to\infty$. Since the map $T$ is not assumed to leave the target invariant, the proposed algorithm is called non-equilibrium orbit importance sampling (NEO-IS).   The authors propose to use conformal Hamiltonian dynamics as a mapping (standard Hamilton dynamics based on the negative log target plus a friction term). The dissipative friction term allows the authors to control how exhaustively phase space is explored.  The second main contribution is an MCMC algorithm based on the ideas used in NEO-IS. NEO-MCMC applies iterative Sampling Importance Resampling (SIR) within the NEO-IS framework. NEO-MCMC jumps from one point on an orbit to a point on another orbit. Orbits are sampled according to their importance weight.  The algorithms are benchmarked on various targets (25 Gaussian mixture, Funnel) and applied to training a VAE.  

The paper puts forward a new family of algorithms, NEO, for computing normalizing constants and sampling from complex distributions. NEO-IS is an importance sampling estimator of the normalizing constant. NEO-MCMC combines NEO-IS unbiased estimator of the normalizing constant with iterated sampling- importance resampling methods.

The paper introduce two sampling approaches named NEO-IS and NEO-MCMC respectively. The former method aims to estimate normalizing constant for an unnormalized target distribution, and the latter method can be executed to do sampling from unnormalized distribution. Empirical experiments reveal impressive performance for both the approaches on Normalizing Constant Estimation and Sampling tasks.  Exact theoretical analysis, especially convergence guarantees are also given in this paper. 

This paper proposes NEO which extends Non-Equilibrium Importance Sampling (NEIS) to use discrete orbit and extra weights. The probability distribution and sampling process is derived based on the Jacobian determinant. This method could be applied to estimate the partition function and a biased estimation called self-normalized IS is proposed based on estimated partition function. This paper also derive upper bound for the bias and MSE. A damped Hamiltonian system is adopted to equip the original system with extra momentum and provide orbit trajectory for NEO. NEO is also combined to Sampling Importance Resampling to derive the NEO-MCMC procedure. The authors present reasonable experiments to evaluate their methods.