BNEM: A Boltzmann Sampler Based on Bootstrapped Noised Energy Matching
Authors: RuiKang OuYang
Bo Qiang
José Miguel Hernández-Lobato
Venue: Accepted by TMLR ()

Abstract: Generating independent samples from a Boltzmann distribution is a highly relevant problem in scientific research, \textit{e.g.} in molecular dynamics, where one has initial access to the underlying energy function but not to samples from the  Boltzmann distribution. We address this problem by learning the energies of the convolution of the Boltzmann distribution with Gaussian noise.  These energies are then used to generate independent samples through a denoising diffusion approach. The resulting method, \textsc{Noised Energy Matching} (NEM), has lower variance and only slightly higher cost than previous related works. We also improve NEM through a novel bootstrapping technique called \textsc{Bootstrap NEM} (BNEM) that further reduces variance while only slightly increasing bias. Experiments on a collection of problems demonstrate that NEM can outperform previous methods while being more robust and that BNEM further improves on NEM. Codes are available at \url{https://github.com/tonyauyeung/BNEM}.

Published in Transactions on Machine Learning Research (03/2026)
A Boltzmann Sampler Based on Bootstrapped Noised Energy Matching
BNEM:
RuiKang OuYang" University of Cambridge
ro352@cam.co uk
Bo Qiang Universily of Washinglon
bqiang.uuw... edu
José Miguel Hernández-Lobato University of Combridge
jmh233 uk
Reviewed on OpenReview: https: // openreveew. net/ forum? id= UOUUU
Abstract
Generating independent samples from Boltzmann distribution is highly relevant problem in scientific research, e.g. in molecular dynamics, where one has initial access the underlying energy function but not to samples from the Boltzmann distribution. We address this problem by learning the energies of the convolution of the Boltzmann distribution with Gaussian noise. These energies are then used to generate independent samples through a denoising diffusion approach. The resulting method, NoISED ENERGY MATCHING (NEM), has lower variance and only slightly higher cost than previous relatod works. We also improve NEM through novel bootstrapping technique called BOOTSTRAP NEM (BNEM) that further reduces variance while only slightly increasing bias. Experiments collection of problems demonstrate that NEM can outperform previous methods while being more robust and that BNEM further improves on NEM. Codes are available at https://githut www//ooyaayouugn///ii...
1 Introduction
A fundamental problem in probabilistic modeling and physical systems simulation is to sample from a target Boltzmann distribution (2) x exp(-E(x)) specified by an energy function E(x). A prominent example is protein folding, which can be formalized as sampling from a Boltzmann distribution (Šledá & Caflisch, 2018) with energies determined by inter-atomic forces (Case et al., 2021). Having access to efficient methods for solving the sampling problem could significantly speed up drug discovery (Zheng et al., 2024) and material design (Komanduri et al., 2000).
Traditional mothods for sampling from unnormalized densitics include Monte Carlo techniques such as AIS (Neal, 2001), HMC (Betancourt, 2018), and SMC (Doucet et al., 2001), which are computionally expensive. Amortised methods using machine learning techniques have been developped in recent years (Noé et al., 2019; Midgley et al., 2023; Vargas et al., 2023a;c; Berner et al., 2024; Albergo & Vanden-Eijnden, 2024; Rissanen et al., 2025). However, most existing methods for sampling from Boltzmann densities have problems scaling to high dimensions and/or are very time-consuming. As an alternative, Akhound-Sadeghh et al. (2024) proposed Iterated Denoising Energy Matching (iDEM), a neural sampler based on denoising diffusson models which is not only computationally tractable but also seems to provide good coverage of modes. Nevertholess, iDEM requires a large number of samples for its Monte Carlo (MC) score estimate to have low variance and large number of integration steps even when sampling from simple distributions. Also, its effectiveness highly depends the choice of noise schedule and score clipping, These disadvantages demand careful hyperparameter tuning and raise issues when working with complicated energies.
"Equllly contributed. Order randomly assigned. Correspondence ro3628ram. and bqiangdue

---

Published in Transactions on Machine Learning Research (03/2026)
Score Reverse SDE: V ( 1)
Update buffer B
time
Outer loop: simulate the diffusion sampler
Inner loop: target energies of noised data
Sample zo
from buffer
time
Optimize
3
Eg
E (2
NEM
Polt
Energy Estimator
Learned Energy
BNEM
D ult
Energy Estimator
Figure 1: Both NEM and BNEM parameterize time-dependent energy network Eo(x,, t) to target the energies of noised data. NEM targets MC energy estimator computed from the target energy function; BNEM targets Bootstrap energy estimator computed from learned energy functions at slightly lower noise level. Contours are the ground truth energies at different noise levels; represents samples used for computing the MC energy estimator, represents samples used for computing the Bootstrap energy estimator, and the white contour line represents the lcarned energy at time u.
To further push the boundary of diffusion-based neural samplers, propose NOISED ENERGY MATCHING (NEM), which learns a series of noised energy functions instead of the corresponding score functions. Despite a need to differentiate the energy network when simulating the diffusion sampler, NEM targets less noisy objectives as compared with iDEM. Additionally, using an energy-based parametrization enables NEM to use bootstrapping techniques for more efficient training and Cedropolis-Hastinne corrections for more accurate simulation. By applying the bootstrapping technique, we propose variant of NEM called BooTSTRAP NEM (BNEM). BNEM estimates high noisc-level encrgies by bootstrapping from current energy estimates at slightly lower noise levels. BNEM increases bias but reduces variance in its training target.
Our methods outperform alternative baselines on both experiments with synthetic and n-body system targets. Additionally, we found that targeting energies instead of scores is more robust, requiring fewer MC samples during training and fewer integration steps during sampling. The latter compensates for the need of our methods to differentiate the energy networks during sampling for score calculation.
Our contributions
follows:
We introduce NEM and prove that targeting noised energics rather than noised scores has advantages through theoretical analysis.
We apply bootstrapping energy estimation to NEM and present thcoretical analysis of the Bias- Vuriance trade-off for the resulting dilfusion-based sampler.
We perform experiments showing that BNEM and NEM significantly outperform alternative baseline neural samplers four different tasks, while also demonstrating robustness under various hyper- parameter settings and efficiency in terms of number of function evaluations.
2 Preliminary
We consider learning a generative model for sampling from the Boltzmann distribution
Ptanget (z) exp(-E(x))
Z = - -()))x,
(1)
Z
2

---

Published in Transactions on Machine Learning Research (03/2026)
where & is the energy function and Z is the intractable partition function. Efficiently sampling from this type of distribution is highly challenging. The recent success of Diffusion Models provides promising way to tackle this problem.
Diffusion Models (DMs) (Sohl-Dickstein et al., 2015; Ho et al., 2020; Song et al., 2020) learn a generative process that starts from a known and tractable base distribution, a.k.a. denoising process, which is the inverse of a tractable noising process that starts from the target distribution. Formally, given samples from the target distribution, 20 Margett: the noising process is an SDE towards known base distribution Pl, which is also known the Diffusion SDE:
dx, = (a tydt + g(t)dur; for E [0,
where f(x,,t t) is called drift coefficient, (t) is the diffusion coefficient and is standard Brownian Motion. Diffusion Models work by approximately solving the following inverse SDE
dr, (f(x,,)) (1) log (x:)) + g())dox,
(3)
where Wn is again standard Brownian Motion and Pt is the marginal density of the diffusion SDE at time t. the example the Variance Exploding (VE) noising process, f(z,,t) = and the perturbation kernel is given by (2)((2)(2))) = (2 Np, xb, of), where of := J g (s) Then the learning objective of DMs is obtained by using Tweedie's formula (Efron, 2011):
LDM (0) E oln (mo (=))))(()
w(t)
Zo
so (z t)
oi
where w(l) is positive weighting function and sa(a,, is score network, parameterized by 0, which targets the conditional scores V Voog Pr)o(za/zo) = V (oogN(x;;; I). Minimising the above objective allows us to approximate the marginal scores V (x)) with (0(
3 Methods
We aim to train a diffusion- siffuusio-bbsse neural sampler that samples from (target(z) = e where we assume we only have access to the energy function & without any known data from the target distribution, as illustrated in Fig. 1. To solve this problem, we introduce our proposed NEM framework and discuss the theoretical advantages of energy matching over score matching. Finally, we describe how bootstrapping (BNEM) can be used to obtain further gains. Full descriptions for training NEM and BNEM are provided in Algs. and 2, and visualized in Fig. 1.
3.1 Denoising iistusion-baaeed Boltzmann sampler
We consider training energy-based diffusion sampler (Salimans & Ho, 2021; Gao et al., 2021) corresponding
to variance exploding (VE) noising process defined by dx, = g(f)dwe, where t E [0, 1], g(t) is a function of time
and w, is Brownian motion. The reverse SDE with Brownian motion w, is dr, = g² Po (x,)d, +g((t
where Pt is the marginal of the diffusion process starting at Po := Mtarget.
Given the energy E(x) and the perturbation kernel qe (x2(2z)) = (Nx of), where of = Sh (s)ds, one can obtain the marginal noised density as
x Po(x) and
Pr(xt) x exp((
= (exp=eeppl
Going a step further, the RHS of Eq. (5) defines a Boltzmann distribution over the noiso- portu-petturbdd distribution Pt. The noised energy is defined the negative logarithm of this unnormalized density
Ex(x) = log EN( (xppp (x)x)) where exp(-&((xx)) x pe(x)).
(6)
For simplicity, we denote V. by throughout the paper.
3

---

Published in Transactions on Machine Learning Research (03/2026)
Training on MC estimated targets We can approximate the gradient of log Pa by fitting score network so( 1) to the gradient of Monte Carlo (MC) estimates of Eq. (6), leading to iDEM (Akhound-Sadegh et al., 2024). The MC score estimator SK and the training objective are
SK(T,, 1): V log
Σ exp( exp(-E((z)1)) K
1
(DEM(TI, 1) (Ie, so (z 1)1
(8)
Alternatively, we can fit an energy network Eg(z,, t) to MC estimates of Eq. (6). The gradient of this energy network w.r.t. input It, Le. (oo(T, t), can then be used to estimate the score required for diffusion-based sampling. The MC energy estimator Eg and the training objective can written as
Ex(x,, == -log
( 1=1
Fofe ~ 21
(9)
K
(NEM(TI, | (Ze, () Eo(xe, 1)1l?
(10)
Notice that SK (x,,,) (ER((T,t)) To enable diffusion based sampling, one is required to differentiate the energy network to obtaín the marginal scores, i.c. V Eo,,,,) t), which doubles the computation of just evaluating Eo(x,, t). Regressing the MC energy estimator Ex does not require to compute the gradient of the target energy & during training but it requires computing the gradient of Eg during sampling. In other words, we change the requirement of having to differentiate the energy function & for training, to having to diffcrentiate the learned neural network Eg for sampling. Our method may have computational advantages when evaluating the gradient of the larget energy is difficult or expensive.
Bi-level iterative training scheme To train the diffusion model on the aforementioned Monte Carlo targets, we should choose the imput locations I, at which the targets are computed. Previous works (Akhound- Sadogh et al., 2024; Midgley et al., 2023) used data points generated by a current learned denoising procedure. We follow their approach and use a bi-level iterative training scheme for noised energy matching. This involves
An outer loop that simulates the diffusion sampling process to generate informative samples These samples are then used to update a replay buffer B.
simulation-free inner loop that matches the noised energies (NEM) or scores (DEM) evaluated at noised versions I, of the samples To stored in the replay buffer.
3.2 Energy-based learning v.s. score-based learning
In this section, aim to answer the following questions.
(Q1) Do NEM and iDEM differ with infinite data, training time, and model capacity?
(Q2) Why might NEM outperform iDEM in practical settings?
To answer Q1, We assume that both NEM and iDEM use fully flexible energy and score networks and that both have converged by training with unlimited clean data. Since both networks are regressing on the MC estimates, at convergence, their network outputs are then the expectations Eg- (xx, ([Ex(x,, and sg- (I,, t) = E[Sx (x,, f)). Noticing that Sx = - -V and by changing the order of differentiation and expectation, we show that NEM and iDEM result in the same score in this optimal learning setting, í.e. 8g- Fgg, in Appendix F.
However, factors such as limited data and model capacity, will always introduce learning error, which is related to the variance targets Zhang et al. (2017); Belkin et al. (2019) and raises Q2. In the following, we characterize Fx by quantifying its bias and variance in order to address this question.
4

---

Published in Transactions on Machine Learning Research (03/2026)
Proposition 3.1. e/exp(-888 exx((=((x)))) (() we have
is Cab-Gaussinn, then with probability 1 over Tofe Na
[ E, 5 exp(E((x))) 20 (xxe log(2/5)
1 K
+0
(11)
VK
where Une (nx = (expemm)nn((xx(((((())]
Proposition 3.1 shows an error bound on the training target used by NEM. This bound is invariant to energy shifts (i.e. adding a constant to the energy) because exp(&((x)) (xx) is invariant under that operation. The bound depends on energy values at I₁ and in the vecinity 2t, i.e. (xy)) and uor()z). In particular, the error bound tends to be small in high-density regions (i.c., low-energy regions) where cxp(E((x;)) is relatively low. This yields reliable energy estimates that may be effectively leveraged for Melropolis-Haslings corrections. The bias and variance of Ex be characterised through the above error bound, provided by the following corollary:
Corollary 3.2. 1f exeppp exp((=(<z())) is sub- Gaussian, then the bias and variance of EK are given by
vor(2e)
(
(aas(Fx((r,,))) 0)] R[Fx (Nee (x), Ec(xt))
+ O
(12)
2mg (2xe) ) K2
(xc( +
1
(ar t)]
(13)
m²(x,)
K2
where mx(x) = exp(ppppppp Ex(()) and vol(x)) is defoned in Proposition 3.1.
The complete proofs of Proposition 3.1 and Corollary 3.2 are given in Appendix C. Additionally, for I-dimensional inputs, the variances of Sx and Ex can be related as follows.
Proposition 3.3. Let x R. If exp(-E((x)) 0, exp(-E(=! sub- sub-Gaussion, and ||V exp(p E(x(-()) is bounded. Then
Vor(Sx (2(,t) 4(1 + (VV Var Ex (I2,
1 K
(14)
Proposition 3.3 shows that the MC energy estimator has less variance than the MC score estimator in the I-dimensional case, resulting in a less noisy training signal. The complete proof for the above results is provided Appendix D. We argue that this variance gap naturally extends to nigher-dimennionall settings, which trivial to show when there is perfoct dependence or independence between data dimensions. Therefore, claim that NEM provides an easier learning problem than iDEM. While NEM uses largets with lower variance, it still requires differentiating the energy Eg during sampling, unlike iDEM, which directly uses the learned score 8g. This makes it hard to quantify the actual differences between these methods theoritically. Despite this, our experiments show that NEM performs better in practice than iDEM, in terms of both performance and convergence rale, indicating that the differentiation of Fg during sampling does not introduce additional significant errors in NEM. To address potential memory issues arising from difforentiating the energy network in high-dimensional tasks, we also introduce memory-efficient sampling procedure in Appendix J.
3.3 Improvement with bootstrapped energy estimation
Using an energy network that directly models the noisy energy landscape has additional advantages. Intuitively, the variances of EK and Sк explode high noise levels as a result of the VE noising process. However, we can reduce the variance of the training target in NEM by using the learned noised energies at just slightly lower noise levels rather than using the target energy at time = Based this, propose Bootstrap
5

---

Let us consider a bootstrapping-oner scenario 50-5577 where we assume that we first learn an optimal energy network s by using the MC energy estimator Eq. (9) with unlimited training data and fully flexible model, resulting in Fgr (2,, s) Ex (1,,,,,))) (I); then the Bootstrapped energy estimator is obtained by Published in Transactions on Machine Learning Research (03/2026)
NEM, or BNEM, which uses novel MC energy estimator at high noise levels that is bootstrapped from the learned energies at slightly lower noise levels. Suppose that Eo(;,)) is an energy network that already provides a relatively accurate estimate of the energy at low noise level we can then construct a bootstrap energy estimator at a higher noise level f > s:
k
(x(,,,t,,;0) -lo Σ 1=1
exp( (x)
o ~ (f
(15)
K
(x,, |EE (2,, t, t; - Eg(x,,) - +)12
(16)
where SG) denotes a stop-gradient operation, preventing backpropagation through 0.
plugging Ee; into Eq. (15). Appendix G.1 shows that the bias of the bootstrapping-once estimator can then characterized as
voc(xx)) + 2m²( Bias[Ex t, s;61) 2m²(rz)
(x4
This illustrates that, theoretically, bootstrapping once from time 8 can quadratically reduce the bias contributed at time > 8, while the bias contributed at time 8 is accumulated. By induction, given chain of bootstrapping 80 ... 8; t, we sequentially optimize the network at (Ia,, 84) by bootstrapping from the previous level. Then the bias at (x,,t) can be charecterized follows:
where 0, Sn = 1 and k. >> 0 is a small constant. suppose Eg is incrementally optimized from t E [0, 1] as follows: if E [8,, 8141), (o t) targets an energy estimator bootstrapped from Vs [st-1, 8t] using Eq. (15). For v0 S i S n and vs E 85], the variance of the boolstrap energy estimator approzimately given by
Proposition 3.4. Given {s:li=o such that o², - S
VarlEx(EE 0)] (sr(
Var Ex
(17)
Vos (Te)
and the bias of EK tx is approzimately given by
Vor (xx)) (x 2m, (xx) Ki+I
Uos; (It) 2m2 i3 (xxx Kj. j=1
(18)
where (xx) = Varsicess.ee--ognneeeed
and my (xx)
for 0O 5y<
A detailed discussion and proof given in Appendix G. Proposition 3.4 demonstrates that the bootstrap energy estimator, which estimates the noised energies by sampling from an me-mean Gaussian with smaller variance, can reduce the variance of the training target, because the variance ratio small, while this new target can introduce accumulated bias.
< is very
Proposition 3.4 shows that, the bias of BNEM consists of two components: (1) the target bias term which is reduced by a factor of K compared NEM, and (2) the sum of intermediate biases, which are each reduced by factors of KJ where j 2 Since K is typically large and i 1, the target bias term in BNEM is substantially smaller than in NEM. While BNEM does introduce additional accumulated terms, these terms are also reduced by powers of K, resulting in an overall lower bias compared to NEM under similar computational budgets. Thercfore, with a proper choice of number of MC samples K and bootstrap trajectory { the bias of BNEM (Eq. (12)) can be smaller than that of NEM (Eq. (18)) while enjoying a lower-variance training target.
Varianco-controlled bootstrap schedule BNEM aims to trade the bias of the learning target to its variance. To ensure that the variance of the Bootstrap energy estimator at & bootstrapped from controlled

---

Published in Transactions on Machine Learning Research (03/2026)
by a predefined bound B, f.e. of -o3 s B, we first split the time range [0,1] with 0 = to < <... <tt =1 such that of.. - o S 8/2; then we uniformly sample and t from adjacent time splits during training for variance control.
Training of BNEM To train BNEM, it is crucial to account for the fact that bootstrap energy estimation at & can only be accurate when the noised energy at is well-learned. Given a partially trained NEM energy network, we introduce a rejection scheme for bootstrap fine-tuning to balance the bootstrapped and original MC estimators. Specifically, at each training step, we compute the normalized losses (NEM of and LNEM (Z,,, s) and define an acceptance ratio proportional to their values. If accepted, we use the boolstrapped estimator; otherwise, we fall back to the original MC estimator. The detailed illustration could be found in Appendix B.2 and Alg. 2.
4 Related works
Flow-based Neural Sampler. Flow AIS Bootstrap (Midgley et al., 2023, FAB) trains Normalizing Flows to minimize the a=2-divergence, i.e. the variance of importance weights, which exhibits mode covering properties. Notably, it employs replay buffer for better exploitation.
Path-measure-based Neural Sampler. One can generate samples by sequence of actions, which defines path of the sampling process. Path-measure-base based methods train neural samplers by minimizing the divergence between the sampling path and its reversal. Inspired by diffusion models (Song & Ermon, 2019; Ho et al., 2020), Zhang & Chen (2022); Vargas et al. (2023a); Richter & Berner (2023) propose to learn the sampling process by matching the time reversal of the diffusion process. Doucet et al. (2022) proposes to learn the time-reversal the AIS path, while Vargas et al. (2023c); Albergo & Vanden-Eijnden (2024) introduce a learnable corrector to make sure the sampler transports between sequence of predofined marginal densities. On the other hand, inspired by reinforcement learning, Bengio et al. (2021) proposes GFlowNetss which amortize MCMC by utilizing local information, such as stuj-)trajectory detailed-balancee, to train generative model. Sendera. et al. (2024) proposes techniques to improve GFlowNets. However, most of these methods require simulation during training, which still poses challenges for scaling up to diiheer-imeiiinnal tasks.
Monte-Carlo-based Neural Sampler. Huang et al. (2023) proposes to estimate the denoising score using a Monte Carlo method. To boost scalability, iDEM (Akhound-
et al., 2024) introduces an off off-policy training loop to regress an MC score estimator. It achieves previous state-of-the-art performance on the tasks below. Notably, this MC estimator is used in Vargas et al. (2023b); Grenioux et al. (2024) and resembles importance wnpirtanca-weighhedd estimator of the Target Score Identity proposed by Bortoli et al. (2024). Woo & Ahn (2024) proposes variant to larget the vector field.
Diffusion-based MCMC Sampler. Chen et al. (2024); Grenioux et al. (2024) combine MCMC with diffusion sampling, by a Gibbs-style sampling procedure which alternates between the clean data space and the noisy data space.
Recent developments. Since the initial submission of this work in early 2025, the field of neural samplers has experienced rapid progress. To account for these recent advancements without disrupting the core narrative of our original study, we provide an extended discussion of contemporaneous literature in Appendix A.
5 Experiments
We evaluate our methods and baseline models on potentials. complete description of all energy functions, metries, and experiment setups is in Appendix L. Supplementary experiments can be found in Appendix M. We provide an anonymous link to our code for implementation reference?.
Datasets. We evaluate all neural samplers 4 different problems: a GMM with 40 modes (d = 2),
paprtiil well (DW-4) potential (d 8), a 13-particle Lennard-Jones (LJ-13) potential (d = 39) and a 55-particle Lennard-Jones (LJ-55) potential (d = 165). For LJ-n potentials, the energy can be extreme
Phttps:/ http:////mmmmmoooo/sssee/ee////e

---

Published in Transactions on Machine Learning Research (03/2026)
Table 1: Performance comparison neural samplers across four target distributions. We evaluate each method using three metrics: data Wasserstein-2 distance (x-W)), energy Wasserstein-2 distance (E-W2), and Total Variation (TV). Each sampler is run with three random seeds; we report the mean + standard deviation. indicates divergent training. denotes cases where only one seed converged, with the best result reported. Best values are in bold.
Target
GMM-40 (d
DW-4
=8)
LJ-13 (d = =39)
LJ-55 (d 165)
Sampler
X-M
E-W
TV
X-W,
E-Ws
TV
X-WV
E-W;
TV
X-W,
E-W;
TV
DDS
15.04 am
305.13 mm
0.96
0.82.am
558.
558.79.m..
0.38
PIS
79.86
0.95
FAB
9.08
47.60arm
0.79
112.70 0.
0.38
IDE
60.49.m
0.
2.80
0.16
0.87
67701
0.06g
1.98
12615.77
0.14
NEM
5.
44.56
0.91
a
1.85
0. gm
0.87
g.
o.03...
2m90
118.58
0.
BNEM
3.662m1
1.87.
2.79
0.49
a.38
0.14
0.86
1.02
o.03...
1.88 230
119.46
0.08
(a) Ground Truth
(b) DDS
(c) PIS
(d) FAB
iDEM
(f) NEM (ours)
(g) BNEM (ours)
Figure Sampled points from samplers applied to CMM-40 poteutials, with the ground truth represented by contour lines. For diffusion-based methods, the reverse SDE integration steps are limited to 100.
when particles are too close to each other, creating problems for estimating noised energies. To overcome this issue, we smooth the Lennard-Jones potential through the cubic spline interpolation, according to Moore et al. (2024).
Baseline. We compare NEM and BNEM to the following recent works: Denoising Diffusion Sampler (DDS, Vargas al., 2023a), Path Integral Sampler (PIS, Zhang & Chen, 2022)), Flow Annealed Bootstrap (FAB, Midgley et al., 2023) and Iterated Denoising Energy Matching (iDEM, Akhound-Sadegh el al., 2024). Due to the high complexity of DDS and PIS training due to their simulation- nature, we limit their integration step when sampling to 100. As shown in prior work (Akhound-Sadegh et al., 2024), all basclines perform well on the GMM benchmark when using 1000 integration steps our methods as well. To increase lask difficulty and better assess robustness, we reduce both the number of integration steps and MC samples to 100 in our evaluation. For other tasks, f.e. DW-4, LJ-13 and LJ-55, we stick to using 1000 steps for reverse SDE integration and 1000 MC samples in the estimators. We notice that the latest codebase of iDEM uses 100 MC samples but employs 10 steps of Langevin dynamics before evaluation in the LJ-55 task. For fair comparison, we disable these Langevin steps and stick to using 1000 MC samples for iDEM. For BNEM, we set B = 0.1 for all tasks. All samplers are trained using an NVIDIA-A100 GPU.
Architecture. We use the same net architecture (MLP for GMM and EGNN for n-body systems, i.e. DW-4, LJ-13, and LJ-55) in DDS, PIS, iDEM, NEM and BNEM. To ensure a similar number of parameters for each sampler, the score network is parameterized by (0 0) = fo(x; 1), the energy network is set to be (x, 1 fo (z, t) + e with a learnable scalar c. Furthermore, this setting ensures 1E(3)-invarian)) for the energy network. Since FAB requires invertible architecture, we use in this method a continuous normalizing flow specified by fa(x, t) so that it has similar number of parameters as the other samplers.
Metrics. We applied the 2- Wasserstein distances on data (x-W2) and energy (E-W)) values as quality metrics. Additionally, the Total Variation (TV) is also computed. TV is computed directly from data in GMM experiments and from interatomic distances in n-body systems. To compute W2 and TV metrics, we use pre-generated ground-truth samples: (a) For GMM, we sample from the ground truth distributionn (b) For DW-4, LJ-13 and LJ-55, we use samples from Klein al. (2023b). Note that we take SE(3) into account
In fact, running any MCMC on top of the generated samples trivial, which can be employed by any of the baselines.
8

---

Published in Transactions on Machine Learning Research (03/2026)
Ground Truth
DDS
FAB
iDEM
NEM (ours)
BNEM (ours)
0.06
0.04
0.03
0.04
0.02
0.1
0.02
0.01
-20
-60
-40
-20
-300
-200
0.8
0.4
0.6
1.0
0.4
0.2
0.5
0.2
0
0
2
6
DW-4
IJ-13
1J-55
Figure 3: Histogram for energy (top) and interatomic distance (bottom) of generated samples.
when calculating x-W2 for n-body systems. For a detailed description of the evaluation metrics, please refer to Appendix L.2.
5.1 Main Results
We report x-W2, E-W2, and TV for all tasks in Table 1. The table demonstrates that by targeting less noisy objectives, NEM outperforms DEM on most metries, particularly for complex tasks such LJ-13 and LJ-55. Fig. 2 visualizes the generated samples from each sampler in the GMM benchmark. When the compute budget constrained by reducing neural network size for DDS, PIS, and FAB, and limiting the number of integration steps and MC samples to 100 for baselines none of them achieve high-quality samples with sufficient mode coverage. In particular, iDEM produces samples that are not concentrated around the modes in this setting. Conversely, NEM generates samples with far fewer outliers and achieves the best performance all metrics. BNEM can further improve top of NEM, generating data that are most similar to the ground truth ones.
For the equivariant tasks, f.e. DW-4, LJ-13, and LJ-55, we compute the energies and interatomic distances for each generated sample. Empirically, we found that smoothing the enorgy-distance function by a cubic spline interpolation (Moore al., 2024) to avoid extreme values (i.e. when the particles are too close) is a key step when working with the Lennard Jones potential. Furthermore, this smoothing technique can be applied in wide range of many-particle systems (Pappu et al., 1998). Therefore, NEM and BNEM can be applied without significant modeling challenges. We also provide an ablation study applying this energy energy-smoothing technique to score-based iDEM. The results in Table 3 suggest that it could help to improve the performance of iDEM but NEM still outperforms it.
Fig. shows the histograms of the energies and interatomic distances on each n-particle system. shows that both NEM and BNEM closely match the ground truth densities, outperforming all other baselines. Morcover, for the most complex task, LJ-55, NEM demonstrates greater stability during training, generating more low-energy samples, unlike iDEM, which is more susceptible to instability and variance from different. random seeds.
8-W2 is susceptible to outliers with high energy, especially in complex tasks like LJ-13 and LJ-55, where an outlier that corresponds to pair of particles that are close to each other can result. in an extremely large value of this metric. We find that for LJ-n tasks, NEM and BNEM tend to generate samples with low energies and result in low E-W,, while iDEM can produce high energy outliers and therefore results in the extremely high values of this metric. We also notice that for LJ-55, the mean value of E-W, of BNEM is similar to NEM, while its results have smaller variance indicating better result. Overall, NEM consistently
9

---

Published in Transactions on Machine Learning Research (03/2026)
107
10a
IDEM
DEM
IDEM-1000 IDEM-100
NEM-1000 NEM-100
BNEM-1000
BNEM-100
105
NEM
BNEM
107
BNEM
X-H1)
2-399p
E-W2
10
10
102-
0.050
10a
0.875
10a
10a
10
0.025
101-
0.850
10:
IDEM
NEM
BNEM
IDEM
NEM
BNEM
IDEM
NEM
BNEM
0
0.25
0.50
0.75
1.00
NEE
NEE
Figure 4: Barplots comparing iDEM, NEM, and BNEM evaluations with 1000 vs. 100 inlegration steps and MC samples on the LJ-13 benchmark.
Figure 5: Number of energy cvaluations vs. E-W2- Left: GMM; Right: LJ-55. The y-axis is on loga- rithmic scale.
outperforms iDEM, while BNEM can further improve over NEM and achieve state-of-the-art performance in all tasks.
Robustness of NEM and BNEM. Table highlights the superior robustness of NEM and BNEM, requiring lewer integration steps and MC samples compared to iDEM and other baselines. We further examine this robustness using the LJ-13 benchmark showcase the advantages of NEM and BNEM. Fig. illustrates how metrics change when reducing integration steps and MC samples from 1000 to 100. A comprehensive comparison across different computational budgels is provided in Table 4 (Appendix M). Results indicate that under constrained budgets, iDEM performance significantly deteriorates in GMM and DW-4 potentials, whereas NEM remains more robust. In LJ-13, both methods experience degradation, yet NEM-100 surpasses 1DEM-100 and even 1DEM-1000 in terms of E-W2 and TV metrics. Moreover, BNEM-100 maintains strong performance, matching iDEM-1000 GMM and DW-4, and notably outperforming 1DFM-1000 the I.J-13 benchmark.
Efficiency of NEM and BNEM. Fig. 5 compares the efficiency of NEM and BNEM against iDEM. On the simple GMM task, NEM reaches the same optimality iDEM but requires 5-10 times fewer energy evaluations. BNEM achieves similar results an order of magnitude faster and delivers better performance. On the more challenging LJ55 task, all three methods require similar energy evaluations for convergence. However, NEM and BNEM yield an order magnitude improvement over iDEM. Notably, BNEM converges more slowly than NEM due to bootstrapping instability but eventually reaches a more stable optimum with significantly lower variance in E-W2. These results align with the variance advantage of energy estimators over score estimators established in Proposition 3.3, suggesting that the lower-variance MC energy estimator offers a more stable training signal. We report the memory and runtime (per-loop) overhead of iDEM, NEM, and BNEM in Tables 7 and 8; see Appendix M.5 for more details.
6 Conclusion
We have proposed NEM and BNEM, two neural samplers targeting Boltzmann distributions. NEM regresses a novel Monte Carlo energy estimator with reduced bias and variance, achieving faster convergence rate and better optimality with more robustness to hyper-paramctors when compared to iDEM. BNEM builds on NEM, employing an energy estimator bootstrapped from lower noise-level data, theoretically trading bias for variance. Empirically, NEM outperforms baselines on 4 different benchmarks, and BNEM improves over NEM
Limitations and future work. Though NEM/BNEM achieve better performance on multiple benchmark sampling tasks, they have to differentiate through a neural network, which can increase memory cost. To mitigate this, we have explored using an alternative estimator reduce the memory trace for energy-basod modeling in Appendices and M.7. Another limitation is that, to stabilize training with extreme energy values, we had to cubic-spline interpolation technique for the LJ-n potentials, which introduce small
10

---

Published in Transactions on Machine Learning Research (03/2026)
biases in the generated samples. Future work could avoid this by using contrastive learning alternative energy training scheme.
Broader Impact
This paper presents work whose goal is to advance the field of Machine Learning. There many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Acknowledgements
We acknowledge Javier Antorán, Wenlin Chen, Louis Grenioux, Jiajun He, Zixing Song, and Yusong Wang (listed alphabetically) for helpful and inspiring discussions. JMHL and RKOY acknowledge the support of Turing Al Fellowship under grant CPP/V02376// This project acknowledges the resources provided by the Cambridge Service for Data- Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3 www.cssd3.oom cam nam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Enginecring and Physical Seiences Research Council (capital grant EP EP/T022159/11 and DiRAC funding from the Science and Technology Facilities Council
References
Rishal Aggarwal, Jacky Chen, Nicholas M. Boffi, and David Ryan Koes. Boltznce: Learning likelihoods for boltzmann generation with stochastic interpolants and noise contrastive estimation, 2025. URL https://arxiv.orge org/ 00846.
Tara Akhound-Sadoogh, Jarrid Rector-Brooks Avishek Jocy Bose, Sarthak Mittal, Pablo Lemos, Cheng- Hao Liu, Marcin Sendera, Siamak Ravanbakhsh, Gauthier Gidel, Yoshua Bengio, Nikolay Malkin, and Alexander Tong. Iterated denoising energy matching for sampling from boltamann densities, 2024. URL
https://arxia: httpas/////vvvooo//////// www.oog///////eooooo/oooooee awwpab/2..2226e222//////
Michael S Albergo and Eric Vanden- Nets: A non-cquilibrium transport sampler. arXiv preprint arXiv:2410.02711 2024.
Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal. Reconciling modern machine-learnigg practice and the classical bias variance trade-off. Proceedings of the National Academy of Sciences, 116(32): 15849-15854, 2019.
Emmanuel Bengio, Moksh Jain, Maksym Korablyov, Doina Precup, and Yoshua Bengio. Flow network based generative models for non-iterative diverse candidate generation. Advonces in Neural Information Processing Systems, 34:27381 273944, 2021.
Julius Berner, Lorenz Richter, and Karen Ullrich. An optimal control perspective on diffusion- basee generative modeling, 2024. URL https://arxiv.org/
01364.
Michael Betancourt. A conceptual introduction to hamiltonian monte carlo, 2018. URL https://arxiv:
prg/abbs/701.. 02434..
Donis Blessing, Julius Bernor, Lorenz Richter, Carles Domingo-Enrich, Yuanqi Du, Arash Vahdat, and Gerhard Neumann. Trust region constrained measure transport in path space for stochastic optimal control and inference, 2025. URL https://arxiv org//////// htg//bbs/2508. 12511.
Denis Blessing, Lorenz Richter, Julius Berner, Fgor Malitskiy, and Gerhard Neumann. Bridge matching sampler: Scalable sampling via generalized fixed-point diffusion matching, 2026. URL https://arxis: org/abs/2603.00m
Valentin De Bortoli, Michael Hutchinson, Peter Wirnsberger, and Arnaud Doucet. Target score matching, 2024. URT. https: : //tpx//xivvvv///////// /ttp///www..paaaa////a/aaaaaaaaa http:/abb2222....oe 08667.
11

---

Published in Transactions on Machine Learning Research (03/2026)
David A Case, H Metin Aktulga, Kellon Belfon, Ido Ben-Shalom, Scott R Brozell, David S Cerutti, Thomas E Cheatham III Vinícius Wilian D Cruzeiro, Tom A Darden, Robert E, Duke, et al. Amber 2021. University of California, San Francisco, 2021.
George Casella and Roger I, Berger. Stalistical Inference, volume 2. Duxbury Pacific Grove, CA, 2002.
Wenlin Chen, Minglian Zhang, Brooks Paige, José Miguel Hernández-Tobalo, and David Barber. Diffusive gibbs sampling. arXiv preprint arXi::202.00 030000 2024.
Arnaud Doucet, Nando De Freitas, and Neil Gordon. An introduction to sequential monte carlo methods. Sequential Monte Carlo methods in practice, pp. 3-14, 2001.
Arnaud Doucet, Will Grathwohl, Alexander G Matthews, and Heiko Strathmann. Score-based diffusion meets annealed importance sampling. Advances in Neural Information Processing Systems, 35:21482-2149% 2022.
Bradley Efron. Tweedie's formula and selection bias. Journal of the American Statistical Association, 106 (496):1602 1614, 2011.
RÂomi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, AurAOlic Boisbumon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, LÂOo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamoniy, levgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer. Pot: Python optimal transport. Journal of Machine Learning Reseorch, 22(78):1-8, 2021. URL 20-451. html. 10ml...html.
org wtggpapers/v22/
Ruiqi Gao, Yang Song, Ben Poole, Ying Nian Wu, and Diederik P. Kingma. Learning energy-basee baaee models by diffusion recovery likelihood, 2021. URL htttss://arat org/abs abs0022 2011200012........
Louis Grenioux, Maxence Noble, Marylou Gabrié, and Alain Oliviero Durmus. Stochastic localization via iterative posterior sampling. arXiv preprint orXiv:2402.0 0e0e2.00755. 10758, 2024.
Aaron Havens, Benjamin Kurt Miller, Bing Yan, Carles Domingo Enrich, Anuroop Sriram, Brandon Wood, Daniel Levine, Bin Hu, Brandon Amos, Brian Karrer, Xiang Fu, Guan-Horng Liu, and Ricky T. Q- Chen. Adjoint sampling: Highly scalable diffusson samplers via adjoint matching, 2025. URL https:
aaaxiv.orgg/bss/5004. 11713.
Jiajun He, José Miguel Hernández-Lobato, Yuanqi Du, and Francisco Vargas. Rue: plug-and-play diffuusion inference-time control and energy-based training, 2026. URL https://arxiv.orr org/ htt///bb/25066 25662650005668.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in newral information processing systems, 33:6840 6851, 2020.
Xunpeng Huang, Hanze Dong, Yifan Hao, Yi-An Ma, and Tong Zhang. Reverse diffusion monte carlo. ar preprint TrXie:2307.0003: 2023.
Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Flucidating the design space of diffusion-based generative models. Advances in neural information processing systems, 35:26565 26577, 2022.
Leon Klein, Andreas Krämer, and Frank Noe. Equivariant flow matching. A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances in Neural Information Processing Systems, volume 36, pp. 59886-59910. Curran Associates, Inc., 2023a. URL https:////ppoooeeednggg neurips cwwrripscccppappfffilees files/ pdf..
Leon Klein, Andreas Krämer, and Frank Noé. Equivariant flow matching, 2023b. URL https://arxiv..rgg abs/2306.1500 15030.
Jonas Köhler, Leon Klein, and Frank Noe. Equivariant flows: Exact likelihood generative learning for symmetric densities. In Hal Daumó III and Aarti Singh (ods.), Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 5361-5370. PMLRR 13 18 Jul 2020. URL https: htttpp://prooceeingg.... mlr. www.r.ppess/vee
htmml.
12

---

Published in Transactions on Machine Learning Research (03/2026)
R. Komanduri, N Chandrasckaran, and LM Raff. Md simulation of nanometric cutting of single crystal of crystal orientation and direction of cutting. Wear, 242(1-2):60- 88, 2000.
Guan-Horng Liu, Jaemoo Choi, Yongxin Chen, Benjamin Kurt Miller, and Ricky T. Q. Chen. Adjoint 22565.
schrödinger bridge sampler, 2025. URL
Laurence Mlinng Midgley, Vincent Stimper, Gregor N. C. Simm, Bernhard Schölkopf, and José Miguel Hernández- Lobato. Flow annealed importance sampling bootstrap, 2023. URL https://arxiv.orr orggabs//////./oe wwwpa2/...... 01893..
J. Harry Moore, Danicl Cole, and Gabor Csanyi. Computing hydration froe energies of small molccules with first principles accuracy, 2024. URL
Radford M Neal. Annealed importance sampling. Statistics and computing, 11:125-139, 2001.
Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu. Boltzmann generators sampling equilibrium states of many-body systems with deep learning, 2019. URL
RuiKang OuYang, Louis Grenioux, and José Miguel Hernández-Lobato. A diffusive classification loss for learning energy-based generative models, 2026. URL /tttps://axxdv..or//aa/////0/
Art B. Owen. Monte Carlo theory, methods and ezamples.
wtt.:/ttwwwens...maaaaaa....e www.andsd/m////////mmmmm.mmm 2013.
Rohit V Pappu, Reece K Hart, and Jay W Ponder. Analysis and application of potential energy smoothing and search methods for global optimization. The Journal of Physical Chemistry B, 102((8):9725 9742, 1998.
Lorenz Richter and Julius Berner. Improved sampling via learned diffusions. arXiv preprint 25X7::2507.0119. 2023.
Severi Rissanen, RuiKang YuYang, Jiajun He, Wenlin Chen, Markus Heinonen, Arno Solin, and José Miguel Hernández-Lobato. Progressive tempering sampler with diffusion. arXiv preprint arXiv:2506.052. 2025.
Tim Salimans and Jonathan Ho. Should EBMs model the energy the score? In Energy Based Models Workshop - ICLR 2021, 2021. URL
Marcin Sendera, Minsu Kim, Sarthak Mittal, Pablo Lemos, Luca Scimeca, Jarrid Rector- Rector-Brooks, Alexandre Adam, Yoshua Bengio, and Nikolay Malkin. Improved ool-po training of diffusion samplers. Advances in Neural Information Processing Systems, 37:81016-8108 81045, 2024.
Pawel Šledž and Amedeo Caffisch. Protein structure- based drug design: from docking to molecular dynamics. Curvent opinion in structural biology, 48:93-102, 2018.
Jascha Sohl-Dickstein, Eric A. Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using noncquilibrium thermodynamics, 2015. URL https://arxiv.org/ arg/abs//503 03585.
Yang Song and Stefano Ermon. Gonerative modcling by estimating gradients of the data distribution. Advances in neurol information processing systems, 32, 2019.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishck Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint ar Xiv:2000 1.101..1150 1314550, 2020.
Francisco Vargas, Will Grathwohl, and Arnaud Doucet. Denoising diffusion samplers, 2023a. URL https:
//arxiv. wttprvv..rr//ba///2.....1 13834. 13834.
Francisco Vargas, Andrius Ovsianas, David Fernandes, Mark Girolami, Neil D Lawrence, and Nikolas Nüsken. Bayesian learning via neural schrödinger-föllee flows. Statistics and Compuling, 33(1):3, 2023b.
Francisco Vargas, Shreyas Padhy, Denis Blessing, and Nikolas Nüsken. Transport meets variational inference: Controlled monte carlo diffusions. arXiv preprint arXiv:2307.01050, 2023c.
13

---

Published in Transactions on Machine Learning Research (03/2026)
Dongyeop Woo and Sungsoo Ahn. Iterated mergy-based flow matching for sampling from boltzmann densities, 2024. URT.
Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals. Understanding deep learning requires rethinking generalization, 2017. URL https://arxiv..rrttt
03530.
Qinsheng Zhang and Yongxin Chen. Path integral sampler: a stochastic control approach for sampling, 2022.
URL https: /ttpas/////xvvooo wttp///////.............. 15141.
Shuxin Zheng, Jiyan He, Chang Liu, Yu Shi, Ziheng Lu, Weitao Feng, Fusong Ju, Jiaxi Wang, Jianwei Zhu, Yaosen Min, et al. Predicting equilibrium distributions for molecular systems with deep learning. Nature Machine Intelligence, pp. 1-10, 2024.
14

---

Published in Transactions on Machine Learning Research (03/2026)
BNEM: A Boltzmann Sampler
Based on Bootstrapped Noised Energy Matching Appendix
A Extended Related Works
In this section, we discuss recent advancements in the fields of neural samplers and energy-based models that have emerged concurrently with this work.
Path-Mcasure-Based Neural Samplers. Recent methods have focused on improving neural samplers using path measures induced by diffusion processes. Havens et al. (2025) proposes Adjoinl Sampling (AS), a scalable improvement over Akhound-Sadegl! et al. (2024) derived from a Stochastic Optimal Control (SOC) perspective. Liu et al. (2025) extends AS to a Schrödinger Bridge setting (ASBS). Furthermore, Blessing et al. (2026) unifies AS and ASBS from a probabilistic standpoint and proposes variant that yields additional performance boosts.
Enhancing Neural Samplers via Annealing. Traditionally, neural samplers have been trained using fixed temperature. However, recent progress has highlighted the benefits of annealing techniques, which begin with an casier-to-sample high-temperature distribution and progressively anneal to a lower-temperature one. Rissanen et al. (2025); Schhund-Sadege et al. (2024) first introduced the integration of anncaling into neural sampler Training. Building on this, Blessing et al. (2025) systematically analyzes the annealing process through a trust-region perspective, where the annealing acts a constraint on the KL divergence between the previous and updated samplers. This approach provides an automated and optimal temperature schodule, eliminating the need for manual luning.
Energy-Based Training. Both NEM and BNEM primarily regress the marginal energy defined by the diffusion process. However, instead of relying solely on the energy estimators (Equations (9) and (15)), one can leverage a replay buffer to jointly optimize data-driven auergy-based losses. Recent works (Aggarwal al., 2025; He al., 2026; OuYang et al., 2026) propose various regularizers for training energy-based diffusion models, which could polentially be adapted to enhance the training of NEM/BNEM, Intuitively, these regularizers encourage greater and more efficient exploitation of the replay buffer.
B Details of Training NEM & BNEM
B.1
Iterated Training for NEM and BNEM
B.2
Inner-loop of BNEM Training
In the training of BNEM, we first use NEM to obtain an initial energy network and then apply training with the bootstrapped estimator. We then aim to use bootstrapping only when the energy network is significantly better at time s than at time t. We quantify this by evaluating the NEM losses at times t and s given a clean data point Zo- In particular, we could compare [x(x.) = (NEM(E,, s) and (xxx)) = LNEM (), where V. (20, o²I) and xt N (Zx. of1). However, a direct comparison of and I, is not reliable because the intrinsic difficulty of the regression problem changes with time: even for an optimal regressor under squared loss, the minimum achievable MSE equals the conditional variance of the regression target, i.c., inf, BIjy 12 E[Var(Y X)]. Along the noising path in Diffusion models, this variance typically
15

---

Published in Transactions on Machine Learning Research (03/2026)
Algorithm Iterated training for (Bootstrapped) Noised Energy Matching
Input: Energy network F Batch size b, Noise schedule of, Replay buffer B, Num. MC samples K, I-step
Diffusion solver Dr
while Outer-Loop do
B Update(B, Dr. Pv b))
while Inner-Loop do
To B,t (0, 1), xy ( oř)
5:
(NE 1) = Ex [z,, 1) - Ex( t)
LENEM (22, EK (2t, t, SG(O)) (,, Ex t)ll?, Ultt, where t, St t 4,+1
7: 0 + Update(0, Vo VoCNEM)
8: end while
end while
increases with the noise level, so a smaller raw loss at time does not necessarily imply that the energy network is better learned at than at t. Heuristically, the conditional variance of the NEM regression target scales with the variance of the noising kernel, which is of. To factor out this time-dependent scale and compare losses roughly "noise-normalized" unit, we instead use (xxx) = (x,, 8)/o? and 44xx)) = LNEM t)/a?. We then apply a rejection-basee based training rule rather than a hard switch smooth this noisy loss ratio and improve training stability, follows.
(a) Given and To, we first noise to and f respectively, and compute the normalized losses defined above, i.e. (x,)) and (xx));
(b) These losses indicate how well the energy network fits the noised energies at different times; we then compute a = min(1, (x..)/(())))
(c) With probability a, we accept targeting an energy estimator at bootstrapped from s and otherwise, we stick to targeting the original MC energy estimator.
Algorithm Inner-loop of Bootstrap Noised Energy Matching training
Input: Network Es, Batch size b, Noise schedule of, Replay buffer B, Num. MC samples K
while Inner-Loop do To B.sample()
t~ U(0 1), It (N of)
4: + arggi:t : [t,, 1+
5: 8~ ((tnnt1,tt))) 2, N(zo, o?)
6: [(x.)) - [Fx(Tr,
EI 8) //////0
7: 4(xx) + E -) Ea (oue,tt]1 1)llt))//////
at
with probability a. (22, t) LENEM
10:
11:
Otherwise,
12:
(Z,, = (NEM(IT, + Up VoLBNEM)
13:
14: end while
16

---

Published in Transactions on Machine Learning Research (03/2026)
C Proof of Proposition 3.1 and Corollary 3.2
(i) Proposition 3.1. If (eep((E((0g)) sub- Gaussian, then with probability 8 over Toje ~ NZ of),
we have
-) Ex )xxo(or)exxpl (xp(e (ogg2)))) +0
(19)
VK
where (xx) = (axees,)tn/(ex(((E(())
Corollary 3.2. If exp( cxp(-E(x11))) is sub-Gaussian, then the bias and variance of EK are given by
voc(xe)
1
[ias[Ex((r, 0)] = B|EK (xEE [)] Ex(21)
+ +0
(20)
(
K3
(xe(
1
[ar m²{x)) +
(21)
K2
where m²(x) = exp(-E((z))) and Vor (x)) is defined in Proposition 3.1.
Proof. Let's define the following variables
my(1)) (Ix) = exp(-E;(x:))
(22)
Ust (xt) = oo
(23)
Notice that Ex is a logarithm of an unbiased estimator and by the sub-Gaussian assumption of e'eep(--E one can derive that Ex is also sub-Gaussian. Furthermore, the mean and variance can be approximated using a Taylor expansion, a lechnique formally known as the Delta Method Casella & Berger (2002). Specifically, if a random variable has mean m and variance v, then its logarithm logX has the following mean and variance:
Ellog X] = log m
+ O Var(log X) O
(24)
2m
m²
Hence, according to the sub-Gaussianess of Ex's negative exponential, have
(x+)
(x = Ex(1)) +
+
(25)
K2
(1c( m2(
(x
+0
(26)
K2
Recall that random variable X is o-sub-Gaussian if
2
P(IX -EX e) < 2exp
(27)
20
for any e >> 0. Equivalently, for any 8 E (0, 1), with probability at least 1 - 6, we have
2 X EX| S v 20 2oglog
(28)
Hence, one can obtain the concentration incquality for Ex by incorporating the nab-Gaussianess. For any 10, with probability 1 5, we have
Vor (2+)
2
Ex (x,, 1) - TE EEEE 1)]l s
log mz (x+)
(29)
17

---

Published in Transactions on Machine Learning Research (03/2026)
By using the above Inequality Eq. (29) and the triangle inequality
[[ - Ex(x)))l S Ex (xt, 1) Fx 1)]l + [ (r+, 1)] Ex(+))
(30)
EK (2,, t) El (re, t)]l + For (zz)
(2, (x,) +0(11/22)
(31)
- 2 Uoe
log 2 +
Une (3++)
+0(1/K3)
(32)
mi (zi K
2x (
200
exp(e, (x4)) 1og(2/8) +0(1/K)
(33)
VK
holds with probability 1 for any 8 (0,1).
C.1 Ex error bound v.s. Sк error bound
To connect the bias of MC energy estimator and the MC score estimator, we introduce the error bound of the MC score estimator Sk, where SK VEK, proposed by Akhound-Sadegl. et al. (2024) as follows
20 log(l
log( (3) + [[
(k(( - S(z,, t)ll
(34)
VK
which also assumes that exp(-E(z(f)) is Gab-Gansiian and further assumes sub- Gaussianess over [[ On the other hand, by the sub-Sausssss assumption of exp(-E(=(f)),, it's easy to show that the constant. term C in Equation 34 is C = 22o (xt).
Therefore, we show that the energy estimator error bound would always be 2(1 + (xx)))) times tighter than the score estimator error bound.
C.2 Invariant of the error bounds
At the end, we show that the above error bounds for both EK and Sx invariant to any energy shift, i.c. a shifted system energy E(x) C with any real value constant C will always result in the same bounds.
The proof is simple. By noticing that
0st (xt))
((x)) -e=20Varrienldc
(35)
(36)
my
(37)
ma(
(38)
the exponential shifted value will be cancelled when calculating Use (xxe (x.))m((z)))
D Proof of Proposition 3.3
Proposition 3.3. Let x R. If exp(-E((x))
0, exp(-E(x!())) sub- Gubs-inssion, and
I/V exp(-E(=1))) is bounded. Then
Vor(Sx (2+,1) 4(1 + (E)))) (k) Var Ex (Iz, t)]
(39)
18

---

Published in Transactions on Machine Learning Research (03/2026)
Proof. Review that Sy can be expressed wmportance-weighted estimator follows:
R Σ V eexp(-E(<x))) R Li=1 exp( -E(x))))
(x(
(40)
According to the assumpsion that || (xx)((=((=)))1] is bounded, let | [Vep(--(<o)) S M, where M >> 0. Since bounded variable sub-Gaussian, this assumption resembles a uub-Gaussiiness assumption of V/ )x
On the other hand, according to the assumption that exp(-&((z)) for any in the support, there exists a constant e such that exp( exp(-E(z)j)) 2 c and thus
K i= V Vexp(-E(<1))
(41)
i= exp( EE(x
Σ V Vexp(-E(<0))))
S M/c
(42)
Ke
Therefore, SK bounded by M/c, suggesting it is sub-Gaussian. Recap that the error bound of Sg given by Eq. (34) follows.
20 (Vd
[K((E, t) ((,,, f)ll
(43)
vK
which suggests that we can approximate the variance of Sx(x,, t), by leveraging its sub- Gaussianess, follows,
Var[Sk (xz, t)]
+0
(44)
m2(x))
K2
Therefore, according to Eq. (26) we derive that
4n (a, 1(1 (m
+
VarlSk (Iz,, 1))) Var[Ex (24, t)]
(x,))
K
①
4(1 + (V+++
(
1 K
(45)
(x.)) + +0
(m.)
Step (i) follows from the asymptotic expansion of the ratio. By factoring out the leading order term 1/K from both the numerator and denominator, the quotient converges to the ratio of the coefficients with a residual error of order O(/KK).
E Bias and Error for Learned Scores
In our diffusion sampling setting, the most essential component is the learned score, which can be given by learning the MC score estimator (Eq. (7); iDEM) directly or differentiating a learned energy network (NEM and BNEM). In this section, we discuss the scores learned by iDEM and NEM, well as their bias and error w.r.t. the ground truth noisy score S.(x;) - (x)). Tn E.1 we first show that in optimal case, where the networks are infinitely capable and we train it for infinitely long, NEM would ultimately be equivalent to iDEM in terms of learned scores. Then we characterize the bias and error w.r.t. Sx(x:):
Bias of Learned Scores (F.2): the optimal scores learned by NEM and iDEM, i.e. B[Ex( Q)}, are biased due to the bias of the learning targets Ex (for NEM) and Sx (for iDEM), It results in inreducible errors to (xx)).
Error of Learned Scores (E.3): the neural network theoretically guaranteed to converge to the optimal value, i.e. ([Ex(x,, t)] for NEM and ( t)] for iDEM, with infinitely flexible network and infinite training time. While in practice, the learned networks are often imprefect. It results in non-zero learning errors in practice, which are theoretically reducible.
19

---

Published in Transactions on Machine Learning Research (03/2026)
E.1 NDEM=NEM in optimal case
Suppose both the score network st(r;; t) and the energy network Eo(ar, are capable enough and we train both networks infinitely long. Then the optimal value they learn at (Iz; t) would be the expected value of their training largets, ie. 86- (1), ( 1) and Eg- (æz, 1) ((Ee())) Notice that
Sk (I,, -v, we rewrite the optimal learned scores given by the energy network as follows:
V Eg. (Iz, = = -V, Bl (Fe, t)]
(46)
(-E(ob7))
=-V,, (1:)
log
exp
z0f of11
(47)
Dit
i=1
K
V (1:
Fale
log
K
Σ exp
-8(x, + <bit)) {
Fop
~(e,,0000031)
@
(48)
K
E (1:K)
Va, log
exp
(E(xx + =
(49)
K
K
1
(-E (z611)
E (1:K) Ools
V log
exp
(50)
K
1=1
B[- V Ex(x,
(51)
Sxx ()
(52)
Therefore, in the optimal case, scores learned by both NEM and iDEM would be equal. Then a question arises as to Why is NEM more favorable theoretically and practically, without concern about the additional computation.". We will address this question after presenting the bias for the learned scores. Notice that the above derivation also shows that, compute the expected value of Sk, one can first compute E[Ex (x2, t)] then differentiate this expectation w.r.l. Tt.
E.2
Bias of Learned Scores
Since the regressing targets are biased, the optimal learned scores biased. To derive the bias term, we assume the sub- Gaussianess of (((()) and leverage Eq. (46). According to Proposition 3.1 and Corollary 3.2, the optimal energy network is obtained as follows
(x4) Eg- (Zx, 1) (xx, + 2m (x:) K
(53)
We first review mx(xx) (Eq. (22)) and vor(xx) (Eq. (23)):
me(xt) exp(-Ei(xx))
(54)
(55)
/ [xpp-22(())N o mx())
(56)
20

---

Published in Transactions on Machine Learning Research (03/2026)
We can derive their derivatives w.r.t. 2a follows:
m,(x,)) = (m(
(57)
vor (xx) VoL / (1)p- (z)(2))
(58)
/ eeppl (enpl22(())V,,NNssoo 2m (maaa
(59)
- 2m, (2) (m)
(60)
-
Vz (xc (m((x))mm(a)))
(61)
(62)
2 / (epp= + 2m (E)
(63)
Therefore, the predicted noisy score al (Ie; given by differentiating Eq. (53) w.r.t. Iz, which is
-Ma, (((()) + Vor (3ee
(64)
2mL (x.))())
1 (1's )m (xx)((x)))))))))) - 20n (1x) (2))
=S,(x;)
(65)
2K
(x))
Var(are) + 2000 (xre) (x0c)
=S,(x))
(66)
2m
(K(K)) (
(x,(
+ (E)(a))
(67)
- I /exp(-2E(x) + +2) ( E((r)
1
1 K
Se(Iz)
(xE(z)) 1dr
(68)
K
/ V - (0 + ())
1 + K
S,(11) + 2K
(69)
Therefore, we derive the bias of learned score as follows:
Va, Ege,, - -)((()))
(70)
1 :/ f
K / (((()) E((s)) ((
K Sa(xe)
28
(71)
Eq. (71) shows that the score bias decreases linearly w.r.t, number of MC samples K. In particular, at low noise level (i.e. small t), the latter term of Eq. (68) or Eq. (69) would vanish as (1) (x(x) is close to E(x) and
(2) N(x,,, of)) learned energy network and the directly learned score network.
close to Ft, resulting in a small bias of the scores given by both differentiating the
For BNEMM, given bootstrapping schedule { with su = 0 and 1, the optimal network is given by
Upe (v+))
Vosy (2re)
Eg. (It, t) = (xt) + )x) Ki+1
(72)
2m3 (3))
21

---

Published in Transactions on Machine Learning Research (03/2026)
where E (s St+1]. Then the learned score is obtained by differentiating Eg--
(x4)
Vasy (xt)
-va, Eg. Eg. (Ir, =V 2x.
E((x))
(73)
2m
(2, Kit
2m3 K3
i+1
+
Σ
1
S((11)
(74)
KJ
j=1
1
exp((22((x) + (x,)) &(,,),,)))))))))))))))) V2E(z)) of
K
exp((=28 + 28, (1))) (V,,e, (2t) v.E(x)) (M of N)x
K
which gives us the bias of learned scores for BNEM.
E.3 Error of Learned Scores
Even though the neural networks can converge to the optimal value in theory by assuming infinite training time and infinitely flexible architecture, they are always imperfect in practice. To measure this imperfectness,, we define the 'learning error" the discrepancy between the actual learned value and the optimal learned value. To address the question Why is NEM more favorable thcoretically and practically, without concern about the additional computation!', learning task against regressing the score estimator SK as the former target is shown to have smaller variance in Proposition 3.3, resulting in smaller learning error.
empirically claim that regressing the energy estimator Ex is easier
Let Ea t) = (x, t)] + (g(,,,)) t) be the learned energy network and (3 t) = (x t)] + es(2, t) be the learned score network, where eg and es are the corresponding learning errors. We also assume that eg is differentiable w.r.t. xy. Then the score given by the energy network is
V Eo (x
t)] Ven,ece = 1)] (e,,)()),
(75)
(76)
Therefore, the errors of the network that affects the diffusion sampling are (x (for NEM) and es(x,, I) (for iDEM). Since learning the energy is easier due to the small variance of its training target, we can assume that (= t) is L-Lipchitz. Therefore, the differentiated learning error of the energy network is controlled by the Lipchitz constant L., i.e., ee, f)ll < T.
In practice, this Lipchitz constant can be small because learning the energy estimator is an easier task compared to learning the score estimator, resulting in L t)) which supports the empirical experiments that NEM performs better than iDEM.
F Ideal Bootstrap Estimator = Recursive Estimator
In this section, we will show that an idcal bootstrap estimator can polynomially reduce both the variance and bias w.r.f. number of MC samples K. Let a local bootstrapping trajectory - I, with Suppose we have access to Eg. Let Eg(x., is learned from targetting estimator bootstrapped from En to approximate &,. Suppose Ez(x,, s) is idcal, i.e. Eg(2,. = EK (æ,, 8, E.). The bootstrap energy estimator (Eq. (15)) at t from s can be written as follows:
Ex(2 8; -log
1
Σ
exp(-EE(xx) 8)),
(77)
KL
1=1
22

---

Published in Transactions on Machine Learning Research (03/2026)
By plugging Eg(x=ft== (8) 8) Ex( 8, u; Eu) into the above equation, we have
(, = - log
expl - Ex (x0(-Ex(xg o/e u; Ew))
(78)
K
i=
K
=loo K
exp( 1ep(--Enp Eu (x(
1u ~ NN(=,, -0;--0))1
(79)
J=1
where (65)
Fuls ~ Nx; (K1) ( and Eu are independent, we can combine these Gaussian-noise injection:
(0) (d? - - T) Notice that
(8) + Vo - o3 10 (9) = x.a = + Vo? - o.. and
Fuls
-06 (5)
K
K
Ex(x,, s; = -log K t=1==1
Σ
exp exp(-Eu(=u(t) - Exa (,;
(80)
Fult NN Fe, (o? - I.) Therefore, the ideal bootstrapping at from over this local trajectory (u - + () is equivalent to using quadratically more samples compared with estimating & directly from u. With initial condition & & and
where ((9)
0, we can simply show that: Given bootstrap trajectory {s:] where 0, 8% 1. For any & E [st, s141] and any 8 E [st-1, with i = O,..., n 1, the ideal bootstrap estimator is equivalent to
Ex 8; (Z,,
(81)
which polynomially reduces the variance and bias of the estimator. In other words, the ideal bootstrap estimator is equivalent to recursively approximating &, with Ex (x., 8) utill the initial condition E.
For simplification in appendix G, we term the ideal bootstrap estimator with bootstrapping n times as the Sequential(n) Estimator, ESea(n), i.c.
Sea(()) (zr, t): Ega+i (I(et)
E
(82)
G Proof of Proposition 3.4
Proposition 3.4. Given bootstrap trajectory {s:) such that o², - o2 _ S K, where 8n 0, 8n and K >> 0 is small constant. Suppose Eg is incrementally optimized from t E [0, as follows: if t E St+1l, Eo t) targets an energy estimator bootstrapped from Vs E [s using Eq. (15). For
vO < i < and Vs E [s]- the variance of the bootstrap energy estimator is approzimately given by
VarlEK t, 0) = Vec(Ie)) Var Ex,,,, t)]
(83)
(xx)
and the bios of (xt,,,,,,;; 0) is approzimately given by
Voe
Voo, (xt)
(x t; 0)] = 2m, (me Kitl
+
(84)
2m2 (xx) K3
where m=(x.)
and (x=) Var weppeeennnnnppppeppee for vo <
Proof. The variance of Ex ( L, 8; 0) can simply derived by leveraging the variance of a sub-Gaubun aub-Gaussian random variable similar to Eq. (26). While the entire proof for bias of Ex (,,,,, organized as follows:
we first show the bias of Bootstrap(1) estimator, which is bootstrapped from the system energy
we then show the bias of Bootstrap(n) estimator, which bootstrapped from lower level noise convolved energy recursively, by induction.
23

---

Published in Transactions on Machine Learning Research (03/2026)
G.1 Bootstrap(1) estimator
The Sequential estimator and Bootstrap(1) estimator are defined by:
K exp( Ex olc. s)),
Ex (z,,t): == log
1
o ~ ~(x (o o3)1)
(85)
K
im
-og k
exp( ( ()))))
Top
(0)
N(x
(86)
1 BB()) (zr, (, log exp(-Eo(es) a/c s)),
# N(z (o -0))
(87)
K
The mean and variance of a Sequential estimator can be derived by considering it as the MC estimator with K² samples:
(xe( and Var(ER Seq(1) (zz, t)) = Vox (xx mēr (2) )²
E B|E
Sea(l) (zz, f)] = Ex(r)) + 2mL (2
(88)
While an optimal network obtained by targeting the original MC energy estimator Eq. (9) at s is
Ubo (3r.)
Eg- (Irs 8) = ElEx (1x,, 8) (ogmm(x))
(89)
2m
Then the optimal Bootstrap(1) estimator be expressed
BB (2i,t,;;0)) - llog
Σ exp (-(-Dog my (2 + K
1
K
(z Uns (Ir alc)
(90)
2m (a
Before linking the Bootstrap estimator and the Sequential one, we provide the following approximation which is useful. Let a, b two random variables and {a,) {b,) corresponding samples. Assume that {bş J are close to 0 and concentrated at ma (hence my closes to 0), while {a are concentrated at ma, then
K
DK xi=l exp( (a + b.))
log K Σ exp(-(a; + b.)) = log K
exp( (x;(-a;)
a:)
(91)
ΣR exp( -a:(-a;)
K 21 exp( (((a--aaa + bb))
log
exp(-ai) + 1og
(92)
Σ/ exp(-a:)
1=1 exp(-a,) (-a,) bi) + ( (b,l?)
log
exp( a.) + log
(93)
K
ΣK exp(-ex
i=1
K
1 11=1 exp(-a;)0g
= loog K
exp( -a,) + log
+ + (6,12)
(94)
EK exp(-aa
1=1
ΣK i== exp( - b;)b;
llog
exp( (xp(-a;))
+ Omax (b,1²)
(95)
K
Di exp( -a;)
iml
log (ap(-a,) my + O [b,13]
(96)
where approximation Eq. (93) applies first order Taylor expansion of e 1 + x around since {b: KK are close 0; while Approximation uses log(1 + x) x x under the same assumption. Notice that when
"We consider minimizing the Ig-norm, f.e. =ar ming a Fa(m,, Ex(mn, Since the target, Ex, is noisy, the optimal outputs are given by the expectation, Fig = F[FK]
24

---

Published in Transactions on Machine Learning Research (03/2026)
(9) Nos (x (1)
is large and Ao = of - oo? I x is small 2m .lt x
j are close to 0 and concentrated at (m
vor(2e)
Therefore, by plugging them into Eq. (96), Eq. (90) can be approximated by
Vo (xt) (xx))
E (1;,,, 20 -log
ms (x a/p
(97)
When is large and of o3 is small, the bias and variance of EK (c 9) are small, then we have
log (=( (8) log 1 K
Σ Sep(-Ex(a))) 8)) EK Seq((2) (Ie, t).
(98)
Therefore, the optimal Bootstrap estimator can be approximated as follows:
E BB( B [1 s;;t = ER Sea(1) (2,, + 2m (m)(
Up (xxe)
(99)
where its mean and variance depend on those of the Sequential estimator (Eq. (88)):
t, 8*) ( Ex(1,) + vor (xx(
Un (x²)
E[Ex B(1) (r,,
(100)
2m (xx)
2m (z, K
vx0( m2 (xx)))
Var B(1) (x,,t,;;0")
(101)
G.2 (ootstrap(1) estimator
Given bootstrap trajectory {s: where so 0 and 8n s, and Eg is well learned at [0, Let the energy network be optimal for S Sn by learning a sequence of Bootstrap(i) energy estimators (i S n). Then the optimal value of (2,, given by B[Ex B(n-1)
([EB(())) (a.. s)]. We are going to show the variance of Bootstrap(n)
estimator by induction. Suppose have:
n
Uos (I
Eg- (Z,, = E.
+
(102)
2m,? (3. K3
1=1
=E(E Seq(n (n (Iis) 8)] +
Vos, (xa)
(103)
2m2 (I, K3
K3
Then for any / E (s, 1], the learning target of Eo(x,, 1) is bootstrapped from 8n = 8,
K
Ex B(()) (xr, = = ----log log
exp(
Eg Eg. (1(() s)),
o N(x;0t, (0 0))1)
(104)
K
n-1 Uns (21
= -loo K
Σ
exp
-E[Esegl
Soq(n 1)
(1) o/6
s)]
(105)
2m (z KJ
j=1
Assume that of - 03 small and K large, then we can apply Approximation Eq. (96) and have
K
Vos, (xxe)
F B(n)
(It, t) = llgg K
Σ
exp
(-E
Ex
(eq(n
1)
(x.
(1)
+
(106)
2m3 (xx) K33
1=1
j=1
25

---

Published in Transactions on Machine Learning Research (03/2026)
The first term in the RHS of the above equation can be further simplified as
log
exp
EEEEEK
(x
(107)
=E
log
1
K
exp
--E
Son(n 1 (x/
s/er sn)
Var Fx
Sen(n
(pp) 8)]
(8)
+o
Var Fi
Son(n-
(2.
1)
(3/
a/c: 8)]
K
2
(108)
K
(1 + =((1))
Ups (2ra)
=12()
log
ΣΣ
Σ exp (-E
(109)
Kn+1
2m3 x2 Kn
where ed l.d. N(O,1) I). Therefore, Eq. (106) can be simplified as
B(m)) E (Iz, 1) = -Ee log
K K
exp
(22 (-8 (222)
+ t()))
Wee, (xc) 2m3 )
(110)
Ka+1
==
Taking the expectation over Eq. (110) yields expectation of the Sequential(n) estimator its RHS, we therefore complete the proof by induction:
E
ER B(n) (IZ. 0)] = (0 (1),
log
K K
exp (-E (Fafe +c(1))
+
0os, (x+)
(111)
Kati
2m (Ix K3
i=1
[
Ex (Sqq(n)) (12,
+ Σ
vus (112)
(112)
2m3 (xte) K2
5=1
(or( = Ex(xx) +
V0s, (x0) 2m3 (xx) Ki
(113)
which suggests that the accumulated bias of Bootstrap(n) estimator given by
(x2, (xx,()] Dbr (xx)
V0s (I:)) 2m (xxe) K
Bias Ex B()()))))))))))())
(114)
(x:)
H Incorporating Symmetry Using NEM
We consider applying NEM and BNEM in physical systems with symmetry constraints like n-body system. We prove that our MC energy estimator Ex is G-invariant under certain conditions, given in the following Proposition.
Proposition H.1. Let the product group SE(3) x Sn + O(3)) and pa be a G-invariant density in IRd. Then the Monte Carlo energy estimator of Eg (2t,t) G-invariant if the sampling distribution zoje N(ooji F2, of) is G-invariant, i.e.,
( g z,, of) = Nog Fz, of)
3This simplication leverages the fact that, X, are i.i.d. bounded nab-Gaussian random variables, then log ΣK D exp( EX,) can be approximated by 'Ellog Eme exp(-x.)) Var((X) O Var(()/R)
26

---

Published in Transactions on Machine Learning Research (03/2026)
Proof. Since po is G-invariant, then & is G-invariant well. Let G acts E IR where g o x gz. (6) Since Fofe (z Ft, ot,of) of) equivalent g o TOfe Na g o Ft, Te, of). Then we have
EK (g o x1,t)
log
1
exp(-E(g o =!!!))
(115)
K
K
=llog K
1
Σ
expl exp(-E(<)))) (x(())))
(116)
(6) ((]()
N(
(117)
Therefore, Eg is invariant to G SE(3) x Sn
Furthermore, Ex (2,,,, s; 0) is obtained by applying a learned energy network, which is G-invariant, to the analogous process and therefore is G-invariant well.
Generalizing NEM
This section generalizes NEM in two ways: we first generalize the SDE setting, by considering a broader family of SDEs applied to sampling from Boltzmann distribution; then we generalize the MC energy estimator by viewing it an importance-weighne estimator.
1.1 NEM for General SDEs
Diffusion models can be gencralized to any SDEs as dr; = f(z1, t)dt + g(t)dur, where 2 and wy is a Brownian motion. Particularly, we consider f(x 1) = -a(1), i.e.
dx, at) + g(t)dww
(118)
Then the marginal of the above SDE can be analytically derived
21 = B(t ( B(t)) J (o(s) B(t) = -- - a(o)ds
(119)
where e ~V(00, I). For example, when g(t) = v B(t) and a(t) (tt), where B(t) monotonic function (e.g. linear) inereasing from 3min to Amax; the above SDE resembles Variance Preserving (VP) process (Song et al., 2020). In DMs, VP can be favor since it constrains the magnitude of noisy data across t; while VE process doesn't, i, and the magnitude of data can explode the noise explodes. Therefore, we aim to discover whether any SDEs rather than VE can be better by generalizing NEM and DEM to general SDEs.
In this work, we provide a solution for general SDEs (Eq. (118)) rather than a VE SDE. For simplification, we exchangeably use (() and B,. Given a SDE as Eq. (118) for any integrable functions a and g, we can first derive its marginal as Eq. (119), which can be expressed as:
0 vt = xo + J. (o()))(s))ddee
(120)
Therefore, by defining y1 B, De we have yo zo and therefore:
y1 = yo + / /.. (o(s))(s))
(121)
27

---

Published in Transactions on Machine Learning Research (03/2026)
which resembles VE SDE with noise schedule (² ((s) We can also derive this by changing variables:
dy: (8) x,dt + 3
(122)
3 (1)a(()) + B- (()(=a(()z)d) + g(t)dw;)
(123)
3 (D))((w))
(124)
which also leads to Eq. (121). Let Pe be the marginal distribution of ye and pt the marginal distribution of St. (6) with Moje ~ Nn, we have
P((v)) exp(-E(w))N
(125)
S((1)) = Vy, (og) V log exp(- ((p(-E()))))) (1))
(126)
(g() llg exp((E((o()) (8)
(127)
Therefore, can learn scores and energies of y/ simply by following DEM and NEM for VE SDEs. Then for sampling, we can simulate the reverse SDE of and eventually, we have To yo-
Instead, we can also learn energies and scores of Ze. By changing the variable, we can have
(r())) = 1 xx) = 3, (c(())
(128)
(z,) = S xx) 8, (((i)
(129)
which provides us the energy and score estimator for xt:
K Σ
(x())) ~ - -log 3r
exp(-E((rj))
(130)
K
K
S1(1:) R - log exp(- EE
(131)
: ~ N(x; 3. v4, (()T)
(132)
Typically, is non-negative function, resulting in B(t) decreasing from 1 and can be close to 0 when t is large. Therefore, the above equations realize that even though both the energies and scores for general SDE can be estimated, the estimators are not reliable at large since 3, can be extremely large; while the SDE of y (Eq. (124)) indicates that this equivalent VE. SDE is scaled by 3,1, resulting that the variance of Ye at large & can be extremely large and requires much more MC samples for a reliable estimator. This issue can be a bottleneck of generalizing DEM, NEM, and BNEM to other SDE settings, therefore developing more reliable estimators for both scores and energies is of interest in future work.
1.2 MC Energy Estimator as an Importance- Seppotance-Weiihted Estimator
As for any SDEs, can convert the modeling task to a VE process by changing variables, we stick to considering NEM with VE process. Remember that the MC energy estimator aims to approximate the noised energy given by Eq. (6), which can be rewritten as:
(z,,) = log / (xp(=E(()))(()1a))) (z of
(133)
log f (xx((-((x)) N (Mn;, Tf
/ 40 (xlp))d
(134)
= exp(-E(x)) N(x of1) (x(
log E Eggi (x|pa))
(135)
28

---

Published in Transactions on Machine Learning Research (03/2026)
The part inside the logarithm of Eq. (135) suggests Importance Sampling technique for approximation, by using a proposal 2oje (2)(2)) Notice that when choosing a proposal symmetric to the perturbation kernal, i.e. (0j{(</<;) = (2; of), Eq. (135) resembles the MC energy estimator we discussed in Sec. 3.2. Therefore, this formulation allows us to develop a better cstimator by carefully selecting the proposal (0/x((/22))
However, Owen (2013) shows that to minimize the variance of the IS estimator, the proposal (oje(r[ should be chosen roughly proportional to (()))ttagee (x), where f(z) exp(-E(x)) in our case. Finding such a proposal is challenging in high-dimension space with multimodal Mtarget A potential remedy can be leveraging Annealed Importance Sampling (AIS; Neal (2001)).
J Memory-Efficient NEM
Differentiating the energy network to get denoising scores in (B)NEM raises additional computations compared with iDEM, which are usually twice the computation of forwarding a neural network and can introduce memory overhead due saving the computational graph. The former issue can be simply solved by reducing the number of integration steps half. To solve the latter memory issue, we propose Memory--ffiiient NEM by revisiting Tweedie's formula (Efron, 2011). Given VE noising process, dz, g()) where w, is the Tweedie's formula can be written as:
Brownian motion and of ==
V Vlogp(())) [[ St
(136)
of f rop(wa )
E [x0
(137)
f xo (x(x))
p(zslzo)po dao
(138)
By revisiting Eq. (5), it's noticable that the noised energy, & (Eq. (6)), shares the same partition function as E, f.e. (x) <x)(=E(=))dz, Vt E [0, 1]. Hence, Eq. (138) can be simplified as follows, which further suggests an MC estimator for the denoising score with no requirement for differentiation
f xo NN To, of exp(-expxxe dry exp(-E,
E E[
(139)
exp (- (E(3roj) - (xt))) zoje
(140)
where Zoje ~ N(x; ofI). Given learned noised energy, can approximate this denoiser estimator follows:
Do(xi,t) exp (- (E(2o)) Eo(x,,)) t)) zoje Do(w,,
(141)
exp (- (Fo(Foo Eo ())) 2oje
(142)
where we can alternatively use Do or Do according to the accessability of clean energy E, the relative computation between ((x) and Fo(x,t), f), and the accuracy of (o(x..).
Experimental results of the memory-effecient NEM could be found in appendix M.7.
K TweeDEM: training denoising mean by Tweedie's formula
In this supplementary work, we propose TWEEDIE DEM (Twe)DEM), by leveraging the Tweedie's formula (Efron, 2011) into DEM, í.e. V, log (xx) = Ep(el [0o - 21) o/o/ TweeDEM similar to the iEFM-VE proposed by Woo & Ahn (2024), which variant of iDEM corresponding another family of generative model, flow matching. However iEFM targets the vector field, while TweeDEM targets the expected value of clean data given noisy data.
29

---

Published in Transactions on Machine Learning Research (03/2026)
Table 2: Comparison between DEM, DDM, and TweeDEM.
Score estimator: wis,, with w; Softmao(())
Sampler], Components
Weight Type
Score Type
DEM(Akhound-Soode et al., 2024)
System Energy
exp(-E(<oft)
System Score
DDM(Karras et al., 2022)
Gaussian Density
(1. Ni 24 of
Gaussian Score
V log (a (1),
(ooopi 2e, of!)
TweeDEM
System Energy
exp( exp(-E(Foj()
Gaussian Score
log Mog.N(2r (= 22, 20,
We first derive an MC denoiser estimator, i.e. the expected clean data given noised data 2₁ at t
/ (rop(zol==)dde
E [2]
(143)
f zo
ge
ge
(00 dro
(144)
(2s)
370 M vo, of exp((=((xo) dao /
(145)
(x+())
where the numerator can be estimated by an MC estimator EN of | [xxx(([((x)) and the denominator can be estimated by another similar MC estimator of lexp( (x))), suggesting we can approximate this denoiser through self-normalized importance sampling as follows
K
exp(- (((()))
Dx(2,,t) :
j=1 exp( E 2ofe
(6)
(146)
(i) Wi
(147)
where Foft ~ NNxx,, T) wy are the importance weights and Dx (zt, t) = E[ro|2y]. Then another MC score estimator can be constructed by plugging the denoiser estimator Dr into Tweedie's formula, which resembles the one proposed by Huang et al. (2023)
(8)
K
(1)
Fofr
5x (x+, ;
(148)
of
where --10-- - resembles the vector fields (x(2z) in Flow Matching. In another perspective, these vector fields can be seen scores of Gaussian, i.c. Iog (x; I),, and therefore Šx is an sepoottanc-weiihtee sum of Gaussian scores while SK can be expressed as an importance-weightee sum of system scores -VE. addition, Karras et al. (2022) demonstrates that in Denoising Diffusion Models, the optimal scores are an importance-weightee sum of Gaussian scores, while these importance weights are given by the corresponding Gaussian density, i.e. (p), = Σ, wi (=fft 21)/a? and w; o (Z of We summarize these three different score estimators in Table 2.
L Experimental Details
L.1 Energy functions
GMM. A Gaussian Mixture density in 2-dimensional space with 40 modes, which is proposed by Midgley et al. (2023). Each mode in this density is evenly weighted, with identical covariances,
40
0
Σ
(149)
40
30

---

Published in Transactions on Machine Learning Research (03/2026)
and the means {eike are uniformly sampled from [-40,40] i.c.
40 Pymm (2)
1
(z; E)
(150)
40
Then its energy is defined by the negative-log likelliooodd i.c.
(GMM(x) log Pymm
(151)
For evaluation, we sample 1000 data from this GMM with TORCH.RANDOM.SEDD. following Midgley et al. (2023); Sahound-Sadegh et al. (2024) as a test set.
DW-4. First. introduced by Köhler et al. (2020), the DW-4 dataset describes system with particles 2-dimensional space, resulting in a task with dimensionality d = 8. The energy of the system is given by the deuble-well potential based on pairwise Euclidean distances of the particles,
&DW (xr) = 2T Σ
a(des - do) + b(dij - do)² + c(dij - do)*
(152)
where a, b, c and dg chosen design parameters of the system, T the dimensionless temperature and dej = || =zyll are Euclidean distance between two particles. Following Akhound-Sadede et al. (2024), we set = 0, b = -4 c 0.9 do = 4 and t 1, and we use validation and test set from the MCMC samples in Klein et al. (2023a) as the "Ground truth" samples for evaluating.
LJ-n. This dataset describes a system consisting n particles in 3-dimensional space, resulting in a task with dimensionality d 3n. Following Akhommd-Sadegh et al. (2024), the energy of the system is given by (xx 113 (x)) (x)c(()) with the Lennard-Jones potential
6
12
ELd
Σ
I'm dey
I' dis
(153)
2T
and the harmonic potential
1 Σ
gose (x)
|0o - xcoM|
(154)
where dij = ||,, æyllı are Euclidean distance between two particles, Tm, T and e are physical constants, 2сOM refers the center of mass of the system and the oscillator scale. We use Tm =1, T= 1, 1 and 0.5 the same Akhound-Sadegh et al. (2024). We test our models in LJ-13 and LJ-55, which correspond to d 65 and d = 165 respectively. And we use the MCMC samples given by Klein et al. (2023a) test set.
L.2 Evaluation Metrics
2- Wasserstein distance W2. Given empirical samples μ from the sampler and ground truth samples v, the
2- Wasserstein distance is defined as:
Wa(u,v) (inf / (I. (( y) y)dddy)
(155)
where x is the transport plan with marginals constrained to and v respectively. Following Akhound-Sadegh et al. (2024), we use the Hungarian algorithm as implemented in the Python optimal transport package (POT) (Flamary et al., 2021) to solve this optimization for discrete samples with the Euclidean distance d(r, v) = | -vl2-- x-W; based on the data and E-W2 is based on the corresponding energy. For n-body systems, i.e. DW-4, LJ-13, and LJ-55, we take SE(3), f.e. rotational and translation equivariance into account. For given n-body point cloud pair (X,Y), instead of using Fuclidean distance, the distance is defined
31

---

Published in Transactions on Machine Learning Research (03/2026)
(X, = mESR)) ||X - (YR + +t)l2, where Kabsch algorithm is applied to find the optimal
rotation and translation.
Total Variation (TV). The total variation measures the dissimilarity between two probability distributions. It quantifies the maximum difference between the probabilities assigned to the same event by two distributions, thereby providing a sense of how distinguishable the distributions are. Given two distribution P and Q, with densities p and g, over the same sample space St the TV distance is defined as
1 TV = 2 //2 1p(x) q(x)]dz
(156)
Following Akhound-Sadohh et al. (2024), for low-dimentional datasets like GMM, we use 200 bins in cach dimension. For larger equivariant datasets, the total variation distance computed over the distribution of the interatomic distances of the particles.
L.3 Experiment Settings
We pin the number of reverse SDE integration steps for iDEM, NEM, BNEM and TweeDEM (see appendix K) as 1000 and the number of MC samples as 1000 in most experiments, except for the ablation studies.
GMM-40. For the basic model fo, we use an MLP with sinusoidal and positional embeddings which has 3 layers of size 128 as well as positional embeddings of size 128. The replay buffer is set to a maximum length of 10000.
During training, the generated data was in the range [-1, 1] so to calculate the energy it was scaled appropriately by unnormalizing by a factor of 50. Baseline models are trained with a geometric noise schedule with min le 5, maxx 1; NEM and BNEM are trained with cosine noise schedule with mmin = 0.001 and max 1. We use K 500 samples for computing the Bootstrap energy estimator ER We clip the norm of Sk, and E, to 70 during training and sampling. The variance controller for BNEM is set to be B = 0.2. All models are trained with a learning rate of 5e - 4.
DW-4. All models use an EGNN with 3 message-passinng layers and a 2-hidden layer MLP of size 128. All models are trained with a geometric noise schedule with minn 1e -, Faax and a learning rate of le for computing SK and Eg. We use K = 500 samples for computing the Bootstrap energy estimator EB We clip the norm of SK, so, and V Eg to 20 during training and sampling. The variance controller for BNEM is set to be 8 == 0.2.
LJ-13. All models use an EGNN with 5 hidden layers and hidden layer size 128. Baseline models are trained with a geometric noise schedule with omin = 0.01 and Fmax = 2: NEM and BNEM are trained with a geometric noise schedule with miin = 0.001 and Tmax = 6.0 to ensure the data well mixed to Gaussian. We use a learning rate of Te 3, K 500 samples for FR; and we clipped SK, 8g and V F, to a max norm of 20 during training and sampling. The variance controller for BNEM is set be 3=0 0.5.
LJ-55. All models use an EGNN with 5 hidden layers and hidden layer size 128. All models are trained with
geometric noise schedule with Gmin = 0.5 and 0max = 4. We use a learning rate of le - 3, K= 500 samples for Eg We clipped Sx and 8e to a max norm of 20 during training and sampling. And we clipped V Eg to a max norm of 1000 during sampling, as our model can capture better scores and therefore a small clipping norm can be harmful for sampling. The variance controller for BNEM is set to be 3 0.4.
For all datasets. We use clipped scores as targets for iDEM and TweeDEM training for all tasks. Meanwhile, also clip scores during sampling in outer-loop of training, when calculating the reverse SDE integral. These settings are shown to be crucial especially when the energy landscape is non-smooth and exists extremely large energies or scores, like LJ-13 and LJ-55. In fact, targeting the clipped scores refers to learning scores of smoothed energies. While we're learning unadjusted energy for NEM and BNEM, the training can be unstable, and therefore we often tend to use a slightly larger Tmin- Also, we smooth the Lennard-Jones potential through the cubic spline interpolation, according to Moore et al. (2024). Besides, we predict per-particle energies for DW-4 and LJ-n datasets, which can provide more information on the energy system. It shows that this setting can significantly stabilize training and boost performance.
32

---

Published in Transactions on Machine Learning Research (03/2026)
Table 3: Ablation Study on applying energy smoothing based on Cubic Spline for iDEM
Energy 7
I.J-55 (d = 165)
Sampler 4
x- x-W24
E-W24
TV4
iDEM
2.077 + 0.0238
169347
2601160
0.165 0.0146
iDEM (cubic spline smoothed)
2.086 + 0.0703
12472 + 8520
0.142 + 0.0095
NEM (ours)
1.898 + 0.0097
118.57 + 106.62
0.0991 + 0.0194
M
Supplementary Experients
M.1 Explaination baseline performance difference LJ55 dataset
In our experiment, we observed a notable difference in performance on the LJ555 dataset compared to the results reported in previous work HAkhound-Sadege et al., 2024). Upon examining the experimental setup, we found that the author the previous study employed an additional round of Langevin dynamic optimization when sampling from the IJ555 target distribution. Specifically, the iDEM method only provides initial samples for subsequent optimization. However, this technique is not utilized with simpler distributions, which could lead to increased compulational burdens when sampling from more complex distributions. And notibly, these Langevin steps could be applied to samples generated by any samplers. To ensure fair comparison, we did not conduct any post-optimization in our experiments.
M.2 Comparing the Robustness of Energy-Matching and Score- Matth-na
In this section, we discussed the robustness of the energy- mnergymmatching model(NEM) with the score-matching (DdM) by analyzing the influence of the numbers of MC samples used for estimators and choice of noise schedule on the sampler's performance.
3.5
1.2
2.2
2.5
3.2
3.5
Log
Figure 6: Comparison of the Energy Wasserstein-2 distance in DW4 benchmark between DEM and NEM across varying numbers of MC samples.
Robustness with limited compute budget. We first complete the robustness discussed in Section 5, by conducting experiments more complex benchmark LJ-13. reports different metrics of each sampler in different settings, i.c. 1000 integration steps and MC samples v.s. 100 integration steps and MC samples, and different tasks.
Robustness v.s. Number of MC samples. As in Fig. 6, NEM consistently outperforms iDEM when more than 100 MC samples are used for the estimator. Besides, NEM shows faster docline when the number of MC samples increases. Therefore, we can conclude that the low variance of Energy matoryyg makes it more beneficial when boost with more MC samples.
33

---

Published in Transactions on Machine Learning Research (03/2026)
Table 4: Neural sampler performance comparison for different energy functions. The number after the sampler, e.g. NEM-100, represents the number of integration steps and MC samples 100. We measured the performance using data Wasserstein-2 distance(x-Wip W2), Energy Wasserstein-2 distance(E-Wg) and Total Variation(TV).
Energy 7
GMM-40 (d 2)
DW-4 (d=
LJ-13 (d 39)
Sampler
X-W
E-W24
TVL
X-W22
E-W24
TV
X-W22
E-W22
TV1
1DEM-1000
4.21
1.63eu.a1
0.81ea.0
0.42
1.80
0.13
0.87
77515.90
10
0.06#0.0z
iDEM-100
25.48
60.49270.1
0.82go.0
0.50
2.80g1.v
0. 06
0.88g0.00
1190.59am0
0.07
NEM-
2.
1.6820.0
0.8 00000
0.46gg.0
0.28
0.28 1a
0.02go.01
5.01 29
0.03 00.00
NEM-100
5. 2.28
44. 54.56
0.9 200
0.4829.0
0.85 52
0. 04
0.88
13. 23.14a 14
0.04
BNEM-1000
20.4
0.36 32
0.66
0.49
0.29 2u.0s
0.86
0.62 sovus
0.03 20.00
BNEM-100
3. 2..00
1.87
0. a..9
0 a0.0
a.38
0. 0.14
0.87 au..ww
5 a.93
0.03 200,00
(a) Geometric: 4.16
(b) Cosine: 6.32
(c) Quadratic: 3.95
(d) Linear: 9.87
(a) Geometric: 2.64
(b) Cosine: 2.29
(c) Quadratic: 2.53
(d) Linear: 3.13
Figure 7: Comparison of different sampler (Above: iDEM; Below: NEM) when employing different noise schedules. The performances of X-W are listed.
Robustness v.s. Different noise schedules. Then, we evaluate the performance differences when applying various noise schedules. The following four schedules were tested in the experiment:
Geometric noise schedule: The noise level decreases geometrically in this schedule. The noise at step given by: o; = db 1-t of where dp = 0.0001 the initial noise level, a = the maximum noise level, and is the time step.
Cosine noise schedule: The noise level follows a cosine function time, represented by: 0t 01 coss where 8 = 0.008 hyper-parameter that controls the decay rate.
Quadratic noise schedule: The noise level follows a quadratic decaysa; = anl where op is the initial noise level. This schedule applies a slow decay initially, followed by a more rapid reduction.
Linear noise schedule: In this case, the noise decreases linearly time, represented a, ai
34

---

Published in Transactions on Machine Learning Research (03/2026)
The experimental results are depicted in Fig. 7. It is pretty obvious that for iDEM the performance varied for different noise schedules. iDEM favors noise schedules that decay more rapidly to O when & approaches O. When applying the linear noise schedules, the samples are a lot more noisy than other schedules. This also proves our theoretical analysis that the variance would make the score network hard to train. On the contrary, all schedules are able lo perform well on NEM. This ralustrates that the reduced variance makes NEM more robust and requires less hyperparameter tuning.
Robustness in terms of Outliers. Based main experimental results, we set the maximum energy (GMM-40: 100; DW-4: 0; LJ--3: O; LJ-55: -150). We remove outliers based these thresholds and rocomputed the E-W2. We report the new values as well as percentage of outliers in Table 5, which shows that the order of performance (BNEM>-EM~>iDEM) still holds in terms of better &-W2 value and less percentage of outliers.
Table 5: 8-W2 w/o outliers (outlier%) for different models and datasets. Bold indicates the best value and underline indicates the second one.
Sampler Energy-
GMM-40 (d = 2)
DW-4 (d = 8)
LJ-13 (d 39)
LJ-55 (d 165)
iDEM
0.138 (0.0%)
8.658 (0.02%)
88.794 (4.353%)
21255 (0.29%)
NEM (ours)
0.069 (0.0%)
4.715 (0.0%)
5.278 (0.119%)
98.206 (0.020%)
BNEM (ours)
0.032 (0.0%)
1.050 (0.0%)
1.241 (0.025%)
11.401 (0.0%)
M.3 Empirical Analysis of the Variance of EK and Sx
To justify the theoretical results for the variance of the MC energy estimator (Eq. (9)) and MC score estimator (Eq. (7)), we first empirically explore 2D GMM. For better visualization, the GMM is set to be evenly weighted by 10 modes located in [-1,1] with identical variance 1/40 for each component, resulting in the following density
VaMM(())
Σ
N
xi Mi:
1
(157)
10
40
while the marginal perturbed distribution at can be analytically derived from Gaussian's property:
10
40 + oo)))
(+(x) = --M- * (1)
Σ 1=1
N
x; p4i
(158)
10
given a VE noising process.
We empirically estimate the variance for each pair of (xt, 1) by simulating 10 times the MC estimators. Besides, we estimate the expected variance over æ for each time t, i.e. Ep,(,))(Va((R) (2,,t)) and (P.(xx)(aar [Var(SK (x,,, t)]]
Fig. 8a shows that, the variance of both MC energy estimator and MC score estimator increase as time increases. In contrast, the variance of Fx can be smaller than that of Sx in most areas, especially when the energies are low (see Fig, 8c), aligning our Proposition 3.3. Fig. 8b shows that in expectation over true data distribution, the variance of Eg is always smaller than that of Sg across E 1].
M.4 Empirical Analysis of the Bias of Bootstrapping
To show the improvement gained by bootstrapping, we deliver an empirical study on the GMM-40 energy in this section. As illustrated in appendix L.3, the modes of GMM-40 are located between [-40, 40j with small variance. Therefore, the sub- Gaussianess assumption natural. According to appendix C and appendix D, the analytical bias of the MC energy estimator (Eq. (9)) and Bootstrapped energy estimator (Fq. (17)) can
35

---

Published in Transactions on Machine Learning Research (03/2026)
t=0.1
tw0.3
t=0.5
tw0.7
tw0.9
tw0.92
tw0.95
t=0.97
two 99
t=1.
(a) Ground truth energy across
[0, 1]
12
Score Estimator (5x)
Variance (x10-53
Energy Estimator (Ex)
0.8
0.4
Dif
0.
07
0.5 Time (t)
I-0.3 1-0.92
t-0.95
t=0.97 t-0.90
L-1.0
(b) Expected variance of estimators
(c) Point-wise variance for & E [0.9, 1]
Figure 8: (a) the ground truth energy of the target GMM from to 1; (b) the estimation of expected variance of from 0 to í = 1, computed by weighted
the variance of estimator at each location with weights equal to the marginal density Psi (c) the variance of MC score estimator and MC energy estimator, and their difference (Var(score]- Varl--aaaaaaaaa far[eneryy]) for from 0.9 to 1, we ignore the plots from t E [0, 0.9] since the variance of both estimators are small. The colormap ranges from blue (low) to red (high), where blues are negative and reds are positive.
be computed by Eq. (25) and Eq. (18), respectively. We provide these two bias terms here for reference,
Unr (Ixz)
Bias( (ix,(Ex(E = 2m, (zz)² )2K (2,,
(159)
voc (2t)
Vos, (zt)
(160)
Kn+1
2m 8; (Ir) K3 K3
j=1
Given a Mixture of Gaussian with K components, po(x) = Ee Th Na; pk, Ex) and E(x) = log me(x) and the can be calculated as follows:
m+(x:) = (xp())
(161)
| N( (xp()xxp--aa (xxpp-l)
(162)
f ofo
(163)
ΣnT TK / ( Ex Ex)dr
(164)
Σ TuN(sei ofl Ex)
(165)
36

---

Published in Transactions on Machine Learning Research (03/2026)
and
Uoz (oe(
(166)
/ N( m mz(
(167)
EnTTH Ty
[x) E,)dæ - -m,(aa)
(168)
j.k
Σ TH / ((IIi2,00 2ig COk (22) 1/2 1/2 dx - mx)
Σja 1/2
(169)
3,k
ES 1/2
Σ
TH Cyh (2T))d/21 2 1/2 22 1/2
N (( Pyk,, o Ezt) - (x)),
(170)
3.k
where
Dig = (2)² + 1ー
(171)
Mik = Ena( (2) Mr + +Ex mk)
(172)
Cys= exp (-2
1 [n En " - PR (E, By + =zT'))) ESn(()) En mal)
(173)
In our GMM-40 case, the covariance for each component are identical and diagonal, i.e. En E = vl. By plugging it into the equations, can simplify the mx(1e) and vor terms follows
m=(xx)= = Σ N(x Mk. (o? + v)
1 K
(174)
Uor (3re)
Σ J.k
1
exp (- 40 (ug (
1 T (MJ ms)) N rri 2 (1s + px), (o? + (2(2Tu)
w/2)r) -m](xs)
(175)
K2
We computed the analytical bias terms and visualize in Fig. 9. Fig. 9a visualizes the bias of the both NEM and BNEM over the entire space. It shows that (1) bootstrapped energy estimator can have less bias (contributed by bias of EK and variance of training target); (2) If Ex is already bias, i.c. the "red" regions in the first row of Fig. 9a when = 0.1, bootstrapping can not gain any improvement, which is reasonable; (3) However, if ER has low bias, i.e. the "blue" regions when = 0.1, the Bootstrapped energy estimator can result. lower bias estimation, superioring MC energy estimator; (4) In low energy region, both MC energy estimator result in accurate estimation. However, in bi-level iterated training fashion, we
and Bootstrapped always probably explore high energy al the beginning. Therefore, due to the less biasedness of Bootstrapped estimator at high energy regions, we're more likely to have more informative pseudo data which can further improve the model iteratively.
On the other hand, we ablate different settings of Num. of MC samples and the variance-control (VC) parameter. We visualize the results in Fig. 9b. The results show that, with proper VC, bootstrapping allows us to reduce the bias with less MC samples, which is desirable in high-dimensional and more complex problems.
M.5 Complexity Analysis
To compare the time complexity between iDEM, NEM and BNEM, we let: (1) In the outer-loop of training, have T' integration steps and batch size B; (2) In the inner-loop, have L epochs and batch size B. Let (B)( be the time complexity of evaluating a neural network w.r.t. B data points, Te Te(B) be the time complexity of evaluating the clean energy w.r.t. B data points, and K be the number of MC samples used. Since differentiating a function ff using the chain rule requires approximately twice the computation as shows that
evaluating f, we summurisc the time complexity of iDEM, NEM, and BNEM in Table principle, in the inner-loop, NEM can be slightly faster than iDEM, while BNEM depends on the relativity between complexity of evaluating the neural network and evaluating the clean energy function.
37

---

Published in Transactions on Machine Learning Research (03/2026)
b-0.22 1H052 b-0.77 1-1.00
VC-0.02
VC-0.00
VC=022
VCT
Blas
30
m/o bootstrap
Num. of MC samples
(a) Bias with VC=0.1
(b) Bias with different VC and Num. MC samples.
Figure 9: Empirical analysis on bias with bootstrapping, on GMM-40.
Table 6: Time complexity of different neural samplers.
Sampler4 Phase->
Inner-loop
Outer-loop (or Sampling)
iDEM
(L(2KTa(B) + NN (B)) O (L(KTx((B) + NN (B))) O
(TT
NEM (ours)
O (2TT
BNEM (ours)
(B) (B)
O (2TTxx(B) (B))
Table 7:
Time comparison (in seconds) of different samplers for both inner-loop (network training) and outer-loop (denoising sampling) across different energies
Energy +
GMM-40 (d 2)
DW-4 (d = 8)
LJ-18 (d = 39)
LJ-55 (d 165)
Sampler
Inner-loop
Outer-loop
Inner-loop
Outer-loop
Inner-loop
Outer-loop
Inner-loop
Outer-loop
IDEM
1.657
1.159
6.783
2.421
21.857
21.994
36.158
47.477
NEM (ours)
1.658
2.252
5.217
8.517
14.646
52.503
17.171
114.601
BNEM (ours)
1.141
2.304
25.552
7.547
68.217
52.396
113.641
115.640
Table 8: Memory comparison (in GiB) of different samplers during sampling across different energies
Sampler 4 Energy
GMM-40 (d (d=2) 2)
DW-4 (d = 8)
LJ-13 (d 39)
LJ-55 (d 165)
iDEM NEM/ DEM/BNEM BNEM (ours)
0.02
0.06
0.08
1.30
0.02
0.11
0.25
4.30
Table 7 reports the time usage per inner-loop and outer-loop. shows that due the need for differentiation, the sampling time, t.e. outer-loop, of BNEM/NEM is approximately twice that of iDEM. In contrast, the inner-loop time of NEM is slightly faster than that of iDEM, matching the theoretical time complexity, and the difference becomes more pronounced for more complex systems such as LJ-13 and LJ-55. For BNEM, the sampling time is comparable to NEM, but the inner-loop time depends on the relative complexity of evaluating the clean energy function versus the neural network, which can be relatively higher.
Table reports the GPU memory overhead of iDEM, NEM, and BNEM, where the network architectures are detailed in appendix L.. The batch size for calculation 1024, 1024, 128, and 128, respectively.
38

---

Published in Transactions on Machine Learning Research (03/2026)
Table 9: Comparison between iDEM, NEM, and BNEM, with similar computational budget.
Energy +
LJ-13 (d = (3)
LJ-55 (d = 165)
Sampler +
x-W2t
E-W24
TV
X-W21
8-W24
TV4
iDEM
0.870
6670
0.0600
2.060
17651
0.160
NEM-500 (ours)
0.870
31.877
0.0377
1.896
11018
0.0955
BNEM-500 (ours)
0.866
2.242
0.0329
1.890
25277
0.113
Table 10: Performance comparison between NEM and ME-NEM.
Energy >>
GMM-40 (d=2) 2)
DW-4 (d= 8)
LJ-13 (d = 39)
Sampler +
X-W24
8-We4
TV4
X-W22
8-W24
TV
X-W22t
8-W2t
TV
NEM (ours)
1.808
0.846
0.838
0.479
2.956
0.14
0.866
27.362
0.0369
ME-NEM (ours)
2.431
0.107
0.813
0.514
1.649
0.164
0.88
20.161
0.0338
M.6 Performance Gain Under the Same Computational Budget
It's noticable that computing the scores by differentiating the energy network outputs, í.e. Va, Eo f), requires twice of computation compared with iDEM which computes so by one neural network evaluation. In this section, we limit the compulational budget during sampling of both NEM and BNEM by reducing their integration steps to half. We conduct experiment on LJ-13 and LJ-55, where we reduce the reverse SDE integration steps in both NEM and BNEM from 1000 to 500. Metrics are reported in Table 9. It shows that with similar computation, NEM and BNEM can still outperform iDEM.
M.7 Experiments for Memory-Efficeet NEM
In this section, we conduct experiments on the proposed Memory-Effiient NEM (ME-NEM). The number of integration steps and MC samples are all set to 1000. ME-NEM is proposed to reduce the memory overhead caused by differentiating the energy network in NEM, which leverages the Tweedie's formula to establish 1-sample MC estimator for the denoising score. In principle, ME-NEM doesn't require neural network differentiation, avoiding saving the computational graph. Though it still requires evaluating the neural network twice (see Eq. (142)), this only requires double memory usage of iDEM and can be computed parallelly, resulting in similar speed of sampling with iDEM. In this section, we simply show a proof-of-concept experiment on ME-NEM, while leaving more detailed experiments as our future work.
Table 10 reports the performance of NEM and ME-NEM, showcasing that ME-NEM can achieve similar results even though leverages another MC estimator during sampling.
M.8 Experiments for TweeDEM
In appendix K, we propose TweeDEM, a variant of DEM by leveraging Tweedie's formula (Efron, 2011), which theoretically links iDEM and iEFM-VE and suggests that we can simply replace the score estimator SK (Eq. (7)) with Sx (Eq. (148)) to reconstruct iEFM-VE. We conduct experiments for this variant with the aforementioned GMM and DW-4 potential functions.
Setting. We follow the ones aforementionedl, but setting the steps for reverse SDE integration 1000, the number of MC samples 500 for GMM and 1000 for DW-4 We set quadratic noise schedule ranging from 0 to 3 for TWe DW-4
To compare the two score estimators Sx and Sx fundamentally, we first conduct experiments using these ground truth estimations for reverse SDE integration, í.e. samplers without learning. In addition, we consider using neural network to approximate these estimators, dee. iDEM and Twe.DEM.
39

---

Published in Transactions on Machine Learning Research (03/2026)
(a) Sx (ours)
(b) Sk
(c) TweeDEM (ours)
(d) DEM
Figure 10: Sampled points from samplers applied to GMM-40 potentials, with the ground truth represented by contour lines. SK and SK represent using these ground truth estimators for reverse SDE integration.
Table 11: Sampler performance comparison for GMM-40 and DW-4 energy function. mcasured the performance using data Wasserstein-2 Distance(x-Wa)) Energy Wasserstein-2 distance(E-Wa)) and Total Variation(TV).. i We compare the optimal number reported by Woo & Ahn (2024) and Akhound-Sadegh et al. (2024). indicates metric non-reported.
Energy +
GMM-40 (d =0=] 2)
DW-4 (d 8)
Sampler
x-W24
E-Wo
TV
X-W2
8-W2
TV
SK
2.864
0.010
0.812
1.841
0.040
0.092
ŠK (ours)
2.506
0,124
0.826
1.835
0.145
0.087
iDEM
3.98
0.81
2.09
0.09
iDEM (rerun)
6.406
46.90
0.859
1.862
0.030
0.093
iFFM-VEi
4.31
2.21
iFFM-OT
4.21
2.07
TweeDEM (ours)
3.182
1.753
0.815
1.910
0.217
0.120
Table 11 reports X-W2, E-W2, and TV for GMM and DW-4 potentials. Table shows that whon using the ground truth estimators for sampling, there's no significant evidence demonstrating the privilege between Sx and SK- However, when training a neural sampler, TweeDEM can significantly outperform iDEM (rerun), iEFM- VE, and iEFM-OT for GMM potential. While for DW4, TweeDEM outperforms iEFM-OT and iEFM-VE in terms of x W2 but are not as good as rerun iDEM.
Fig. 10 visualizes the generated samples from ground truth samplers, i.e. Sx and Šx, and neural samplers, d.e. TweeDEM and iDEM. It shows that the ground truth samplers can generate well mode-concentrated samples, as well as TWe,DEM, while samples generated by iDEM are not concentrated on the modes and therefore result the high value of W2 based metrics. Also, this phenomenon aligns with the one reported by Woo & Ahn (2024), where the iEFM-OT and iEFM-VE can generate samples more concentratod on the modes than iDEM.
Above all, simply replacing the score estimator SK with Sx improve generated data quality and outperform iEFM in GMM and DW-4 potentials. Though TweeDEM can outperform the previous statc-of-the-art sampler iDEM on GMM, it is still not as capable as iDEM on DW-4 Except scaling up and conducting experiments on larger datasets like LJ-13, combing Sx and Šx is interest in the future, which balances the system scores and Gaussian ones and can possibly provide more useful and less noisy training signals. In addition, we are considering implementing a denoiser network for TweeDEM as our future work, which might stabilize the training process.
40

---

Published in Transactions on Machine Learning Research (03/2026)
Table 12: Experimental results for GMM and DW-4 compared with Hamiltonian Monte Carlo (HMC) and Parallel Tempering (PT). The number of evaluations of the energy function is indicated in parenthescs. HMC and PT evaluate the true target energy, while NEM and BNEM use learned energy functions.
Energy +
GMM-40 (d=2)
(d=2)
DW-4 (d (d=
Sampler 4
x-W24
E-W2t
TV
x-W24
E-W21
TV
HMC(100)
12.40
5.73
0.87
0.56
455.00
0.45
HM0()100)
11.26
0.03
0.85
0.42
3.44
0.11
PT(100)
9.02
1945.67
0.79
0.47
38.53
0.20
PT(1000)
5.24
2.17
0.85
0.49
3.70
0.10
NEM((10)
5.28
44.56
0.91
0.48
0.85
0.14
BNEM(110)
3.66
1.87
0.79
0.49
0.38
0.14
M.9 Experimental results on non-neural samplers
As shown in Table 12, we compared our methods with traditional samplers that necessitate direct evaluation of the larget energy function during sampling. We conducted hyperparameter sweep the step sizes for both methods based E-W2. The results presented in Table 12 represent the mean performance from three runs with varying random seeds.
In the GMM-40 experiments, Hamiltonian Monte Carlo (HMC) demonstrates impressive performance in sampling low-energy data points. However, it struggles to accurately weigh the different modes that yield high x-W2. Merely increasing the number of sampling steps does not address this issue, indicating that HMC often becomes trapped in local minima when dealing with sharply defined target distributions, On the other hand, introducing more randomness into the sampling process significantly increases the number of outliers. For Parallel Tempering (PT), incrcasing the number of sampling steps results in lower energy values and more accurale capture of the modes of the target distribution. However, even with ten times the sampling steps of BNEM, it still underperforms compared to our best neural sampler. In the DW-4 experiments, both PT and HMC demonstrate good performance on x- W22 however, they still produce a significant number of outliers, which contribule high energy.
41

---

Published in Transactions on Machine Learning Research (03/2026)
Assets and Licenses
Datasets
DW-4, IJ-13, and IJJ55 (Klein et al., 2023b): MIT license
Codebases
DW-4, IJJ13, and LJ-55 target energy functions: MTT license
FAB (Midgley et al, 2023): MIT license (tttprr/ggthab.ccmm/lllla/laabbb
iDEM (Akhound-Sadegh et al., 2024): MIT license (ftbps///ghhuub.ccm//araa////)))
42


Reviews:
[{"review_id": "3RsaXZ5uZt", "reviewer": "Reviewer_Rcqo", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes NEM that reduces the variance and BNEM which uses a bootstrapping target to further reduce the variance at the cost of increasing the bias. Both approaches learn the energy function as opposed to iDEM that learns the score function instead. A strength is that the approach indeed works with reduced variance that is quite strongly validated by experiments. However, it is not supported well by the theory and math that come with confusing notations and the following questions.", "claims_and_evidence_yn": "No", "claims_explanation": "As explained above, it seems that the claims are not supported well by the theory for which I addressed the issues below. In the review, I focused on the technical aspects especially the math in the body and the proofs.", "audience_yn": "Yes", "audience_explanation": "It is based on the diffusion methods that are actively studied in ML; the Boltzman sampler would be an interesting topic in probabilistic ML.", "requested_changes": "# Main Concerns\n\n* **Confusing Notations.** The authors introduced in Eqn. (5) the Gaussian distribution/pdf $\\mathcal{N}(x; x_t, \\sigma_t I) \\equiv \\mathcal{N}(x_t; x_0, \\sigma_t I)$ and used it thereafter. I understand that they are equivalent as $x = x_0$ implicitly and $\\|x_t - x_0\\|^2 = \\|x_0 - x_t\\|^2$. However, it is still confusing in the subsequent formula. For example, \n   * $p_t(x_t)$ in Eqn. (5) is given w.r.t. $x_t$ but $\\mathcal{E}_t(x)$ in Eqn. (6) is given w.r.t. $x$ which is the integration variable in Eqn. (5). The expression (6) is rewritten as $\\mathcal{E}_t(x) = - \\log p_t(x)$, but what is the difference between $p_t(x)$ and $p_t(x_t)$ according to Eqn. (5)? Eqn. (5) depends on both $x_t$ and $x$.\n   * $x_{0|t}^{(i)} \\sim \\mathcal{N}(x; x_t, \\sigma_t^2 I)$ in Eqns. (7), (9), etc., but what does it mean by $x$? It'd be better understandable if it is written as $x_{0|t}^{(i)} \\sim \\mathcal{N}(x_{0|t}^{(i)}; x_t, \\sigma_t^2 I)$ or simply $x_{0|t}^{(i)} \\sim \\mathcal{N}(x_t, \\sigma_t^2 I)$.\n   * Moreover, the notation $x_{s|t}^{(i)} \\sim \\mathcal{N}(x; x_s, (\\sigma_t^2 - \\sigma_s^2) I)$ in Eqns. (15), etc. is inconsistent unless $x_s$ is replaced by $x_t$. Is $x_s$ an errata?\n\n   I believe that those issues arise as the authors always used $x_t$ instead of $x$. I understand that $x$ is a general variable but $x_t$ is a random variable generated by the SDE. E.g., Eqn. (5) does not need to depend on $x_t$ which the authors can replace to $x$ (if $x$ is not used in $\\mathcal{N}(\\cdot)$) and say that $p_t(x)$ is the (marginal?) distribution of $x_t$. Another issue with Eqn. (5) is that $\\mathcal{N}(x; x_0, \\sigma_t^2 I)$ becomes a Dirac delta when $t \\to 0$ (as $\\sigma_0 = 0$), but what is $\\mathcal{N}(x_t; x_0, \\sigma_t^2 I)$ when $t \\to 0$? (both mean and variance converge to zero in this case).\n\n\n* **Technical Points.**\n   - Right above Eqn. (2), the authors said that the base distribution $p_1$ is known, but $p_t(x)$ (at $t = 1$) is a complex integration (5), e.g., not Gaussian. What does it mean by a *known* distribution $p_1$? In general, $p_0$ is not necessarily Gaussian, so that $p_t$ would be Gaussian but only in the limit $t \\to \\infty$, i.e., $p_1$ can be a complex distribution that is hard to be known.\n\n   - $x_{s|t}^{(i)} \\sim \\mathcal{N}(x; x_s, (\\sigma_t^2 - \\sigma_s^2) I)$. Although the variance \"$\\sigma_t^2 - \\sigma_s^2$\" is smaller, the mean $x_s$ is still a random variable which indeed further increases the variance of $x_{s|t}^{(i)}$. Have the authors properly addressed it? If not, then the mathematical claim of BNEM would be just a result of not dealing with that missing variance, although the experiments show that BNEM would reduce the variance.\n\n   - **Bi-level iterative scheme on pp. 4 and Fig. 1.** Although Akhound-Sadegh et al. (2024) and Midgley et al. (2023) used this technique, it is a potential issue that the replay buffer will have a different distribution if $\\nabla E_\\theta$ is not trained well/enough? How do you ensure that the distribution of the replay buffer $\\mathcal{B}$ is the same as $\\mu_\\mathrm{target}$? (I haven’t followed those references though)…\n\n   - **Approximations.** The authors have used equalities \"$=$\" instead of approximations \"$\\approx$\". E.g. Eqns. (17) and (18) in Proposition 3.4 and the equality for the Bias[ . ] on pp. 6 of bootstrapping-once scenario are actually all approximations, according to their proof in the Appendix; the bound (11) in Proposition 3.1 comes with $\\mathcal{O}(1/K)$ according to its proof. The same goes in the proofs (i.e., using $\\approx$ instead of $=$; Eqns (28), (29), (42) need \"$+ \\mathcal{O}(1/K^2)$\"). The authors need to carefully address this for rigorousness and further discuss those approximation preferably in the appropriate place of the main body.\n\n   - How did the authors obtain the Taylor expansions (24) and (25)? Can you give some more details?\n\n   - When Eq. (44) is transformed to Eq. (45), where is the variance $\\sigma_t$ in Eq. (45)? In (41)--(49), is $x_t$, which is also an argument of $E_{\\theta^*}(\\cdot)$, already a random variable or just a general variable (like $x$)?\n\n   - Why the bound in Proposition 3.1 is true with probability $1 - \\delta$? I can't find any relevant proof in the Appendix. \n\n   - When you apply techniques relevant to sub-Gaussians, can you give some more explanations in Appendix, with a proper references? E.g., right above Eq. (24), why $E_K$ is also sub-Gaussian? Also, for Eq. (26). Right above (40), \"Therefore, $S_K$ is bounded by $M/c$, suggesting it is sub-Gaussian.\" --- how?\n\n   - Right above Eq. (95), \"When $K$ is large and $\\sigma_s^2$ is small...\": Is it $\\sigma_t^2 - \\sigma_s^2$? If not, how is $\\sigma_s^2$ small when $s$ is close to 1?\n\n   - How the approximation (95) is valid under the assumption, according to the definitions of $m_s$ and $E_K$?\n\n# Other Comments\n* When Propositions, Lemmas and Corollaries are referred, they are always named \"Theorems\". E.g., right below Proposition 3.1, it is said \"Theorem 3.1 shows an error bound ...\".\n* The numbers of the Propositions and Corollaries and Equations in the Appendices are not matched with the original numbers. E.g., on page 15, Proposition 3.1 and Eq. (11) become Proposition B.1 and Eq. (19).\n* In \" ... the bias of BNEM (Eq. (20)) can be smaller than that of NEM (Eq. (18)) ... \" on pp. 6, \"Eq. (20)\" seems an errata .\n* I understood it, but it'd be good if the authors add more details about why the bound $\\sigma_{t_{i+1}}^2 - \\sigma_{t_i^2} \\leq \\beta/2$ is chosen on pp. 6, rather than $\\sigma_{t_{i+1}}^2 - \\sigma_{t_i^2} \\leq \\beta$.\n* Can the authors explain why Proposition 3.1 is restricted to 1-dim. on pp. 5? What does it mean here by the dimension, and how is it different between 1-dim and n-dim?\n* In the first line on pp. 5, an errata: vecinity $\\to$ vicinity.\n* Can we get any general insight from the results? The Boltzmann distribution (1) can be considered as just another expression of an arbitrary pdf since $\\mathcal{E}$ is arbitrary. Or, is there any restriction on $\\mathcal{E}$ for it to be an energy function?\n* Right above Eq. (88), \"$b_i$ are close to $0$ and concentrated at $m_b$. I believe it also means $m_b$ is close to $0$. If it does, please add one sentence for it.", "broader_impact_concerns": "N / A.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "0fPBdGbn0A", "reviewer": "Reviewer_d83M", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes Noised Energy Matching (NEM) and a bootstrapped variant, i.e., bootstrapped NEM (BNEM), that aim to generate independent samples from a Boltzmann distribution.  Specifically, they work in a setting with access only to an energy function, without samples, and propose an approach that utilizes diffusion-based neural samplers.  Compared with other methods, e.g., Iterated Denoising Energy Matching (iDEM), this work aims to learn noised energy functions rather than score functions.  The authors show that in ideal cases (infinite data, time, and model capacity), NEM achieves the same score as iDEM.  In addition, in practical settings, they theoretically justify reduced variance over iDEM.   Empirically, the proposed methods consistently outperform existing baselines.  \n\nI am not an expert in diffusion, generative models, or energy-based methods.  However, I found the paper to be reasonably clear and accessible. The authors do a good job motivating the problem, explaining the high-level ideas behind their methods, and structuring the theoretical and empirical results in a way that a non-specialist can follow.  For that reason, I'd say the paper is generally well written and clear.  Beyond this, I am not familiar enough with the area to provide a substantive assessment of its strengths and weaknesses.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The empirical claims in the paper are supported by clear numerical experiments comparing the proposed methods with relevant baselines, which are consistent with the authors’ stated contributions. My evaluation is based primarily on the empirical and methodological presentation, as I am not a domain expert and therefore did not evaluate the correctness of the detailed theory.", "audience_yn": "Yes", "audience_explanation": "Generative models are of significant interest to a large portion of the ML community, so this paper would easily be of interest to the TMLR community.", "requested_changes": "I have no specific revision requests.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "NzwyD1wjd8", "reviewer": "Reviewer_6wXU", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces Noised Energy Matching (NEM), a diffusion‐based sampler that learns noised energy functions rather than noised scores, and Bootstrap NEM (BNEM), which further reduces variance by bootstrapping high-noise energy estimates from lower-noise estimates. The authors present theoretical analyses demonstrating lower variance for energy-based Monte Carlo estimators, provide a bias–variance characterization of the bootstrap mechanism, and show strong empirical results on several Boltzmann sampling tasks, including Lennard–Jones systems with up to 55 particles.", "claims_and_evidence_yn": "No", "claims_explanation": "Mostly yes, with a few areas where additional evidence would strengthen the work.\n\n1. Bias introduced by BNEM is not directly quantified. The theory shows accumulated bias, but the experiments only evaluate final sampling quality, not the bias of the learned energies. This weakens the empirical support for claims about bias–variance tradeoffs.\n\n2. Memory and computational overhead of differentiating the energy network during sampling is not evaluated. Since the paper explicitly states this as a limitation, evidence comparing memory or runtime relative to iDEM would strengthen the claims of efficiency.\n\n3. The rejection scheme for bootstrapping is heuristic and lacks empirical justification. Although BNEM works well in practice, ablations or sensitivity studies would make the claim of its effectiveness more convincing.", "audience_yn": "Yes", "audience_explanation": "NEM and BNEM introduce conceptually simple yet impactful improvements to diffusion-based Boltzmann samplers, supported by both theory and experiments. Researchers working on diffusion models, statistical mechanics, and sampling algorithms will likely find the insights especially relevant.", "requested_changes": "1. Report memory and computational overhead during sampling.\n2. A short discussion of the theoretical motivation of the rejection scheme", "broader_impact_concerns": "No", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)