W1-1: No analysis on computational complexity.

To answer your question, we present a table of computational complexity, which includes the number of energy evaluations required, the wall-clock time per training epoch, and the memory footprint on the discrete Ising model 10 x 10, $\beta =0.4$.

W1-2: The bias and the variance of self-normalized importance sampling are not formally discussed.

We provide a formal analysis addressing both the bias from self-normalized importance sampling (SNIS) and the surrogate energy approximation.

Specifically, we show that (1) the bias and variance of the SNIS estimator decay as $O(1/K)$, and (2) with a well-trained surrogate energy, the bootstrapping estimator closely matches the bias and variance of the bootstrapping estimator using the true energy. Moreover, bootstrapping reduces the variance compared to the non-bootstrapped estimator.

We state the propositions with proof sketches; for the second proof, see our response to Reviewer meoi (W1). A full rigorous proof will be included in the revised version.

Proposition. Consider an unnormalized density $\tilde p_1(x_{1})$ and a conditional generator $F_{t|1}^{x_1}(x_{t})$ evaluated on a sample $x_1 \sim q_{1|t}$. Suppose $\tilde p_1(x_1)$ and $\Big\lVert \tilde p_1(x_1) F_{t|1}^{x_1}(x_t) \Big\rVert$ are sub-Gaussian. Then there exists a constant $c(x_t)$ such that with probability at least $1-\delta$,

where $\hat F_t(x_t)$ denotes the SNIS estimator from eq. 9 using $K$ samples $x_1^{(1)}, \\ldots, x_1^{(K)}\sim q_{1|t}$:

Moreover, its bias and variance are expressed as:

proof sketch) By Hoeffding's inequality for sub-Gaussian random variables, there exists a constant $C$ such that

with $1 - \delta$ probability. (Here, $C = \sqrt{\text{Var}[\tilde p_1(x_1)]}$ is a possible choice.)

Since we have bounds for both numerator and denominator of $\hat F_t$, we can also bound the error of $\hat F_t$ to $F_t$. The result is as follows:

Now, let $\hat{A}$ and $\hat{B}$ denote the numerator and the denominator of $\hat F_t$, respectively. We also let $\mathbb{E}[\hat{A}] = A$ and $\mathbb{E}[\hat{B}] = B$. For sufficiently large $K$, the Taylor expansion is expressed as:

To derive the final equation for the bias term, one can express the bias term as follows:

Since $\text{Cov}[\hat{A}, \hat{B}] = \text{Cov}\left[\tilde p_1 F_{t|1}, \tilde p_1\right] / K$ and $\text{Var}[\hat{B}]=\text{Var}[\tilde p_1]/K$, one obtains the conclusion for $\text{Bias}[\hat F_t]$.

To derive the final equation for the variance term, one can combine the sub-Gaussianity of $\hat F_t$ with the error bound on $\left\lVert \hat F_t(x_t) - F_t(x_t) \right\rVert$.

Proposition. Consider an fully trained surrogate energy model $\mathcal{E}_r^{\phi}(x_r) = \mathbb{E}[\hat{\mathcal{E}}_r(x_r)]$, where $\hat{\mathcal{E}}_r(x_r)$ is a biased energy estimator (from eq. 20). Let $\hat F_t(x_t; \mathcal{E}_r)$ be the bootstrapping estimator using $\mathcal{E}_r$ as the intermediate energy function:

Then, the bias and the variance bound of the bootstrapping estimator with surrogate model $\hat F_t(x_t; \mathcal{E}_r^\phi)$ are given by:

where $\hat F_t(x_t)$ is the estimator without bootstrapping. Since $\text{Var}[\tilde p_r(x_r)] < \text{Var}[\tilde p_1(x_1)]$, the variance of bootstrapping estimator is smaller than the $\hat F_t(x_t)$.

W1-4, Q3: Performance on high-dimensional real-world problems is unknown.

While the current experiments are in relatively low-dimensional settings, the Ising model is a widely used benchmark with relevance to real-world physics problems. More broadly, energy-driven training is still at an early stage, and much of the existing literature focuses on validating methods on synthetic or controlled tasks. Scaling to high-dimensional real-world problems is a promising direction for future work in this field.

W2-1: A bit complicated notation $\tilde{p}, \tilde{w}$ and multiple kernels appear in rapid succession, which may force readers to backtrack.

Thank you for the helpful comment. We understand that the use of notations such as $\tilde{p}$ and $\tilde{w}$, along with multiple kernels introduced in close succession, can hinder readability. In the revised version, we will make an effort to briefly clarify these notations whenever they reappear, so that the exposition remains accessible without requiring readers to backtrack.

W2-2: Some important results (intractability of the backward kernels) are primarily deferred to the appendices.

Thank you for the feedback. We agree that the intractability of the backward kernels is important, and we will move the key parts of this analysis from the appendix to the main text in the revised version to improve clarity.

W3-1, Q2: No comparisons with other cutting-edge baselines (e.g., any diffusion-based schemes or generator matching)

We agree on the importance of more baselines; in response, we include results on the continuous many-well 32 benchmark comparing our method with diffusion-based samplers iDEM and PIS:

Continuous task (manywell 32):

Regarding Generator Matching (GM) and Flow Matching (FM), they typically assume access to samples from the target distribution. As such, they are not directly applicable to our problem formulation. We believe that the included baselines (iDEM and PIS), which also operate in energy-only settings, offer a more relevant and fair comparison and help address your concern.

W3-2: Some may view EGM as an incremental extension of GM.

While it is true that our method builds on known components such as importance sampling and generator matching, we believe that their combination in the energy-only setting with bootstrapped estimation represents a significant and broadly applicable contribution.

Unlike generator matching that requires access to data, EGM operates entirely from energy evaluations, making it applicable to a wide range of domains where sampling is difficult or impossible — including discrete, continuous, and hybrid state spaces. In that sense, we view EGM not as a narrow specialization, but as a general framework for neural sampler training across arbitrary state spaces and stochastic processes (diffusions, flows, and jumps).

W4-1: The discussion of connection to other methods (GFlowNet) is presented, but somewhat limited (e.g., no numerical comparison).

Thank you for the suggestion. While our main focus is on proposing a continuous-time neural sampler for general state spaces, we agree that a deeper investigation into the connection with GFlowNets would enrich the broader context of our work. We consider this an interesting and promising direction, and we would like to explore both theoretical and numerical comparisons with GFlowNet-based methods in future work.

Q1: Method's efficacy depending on the choice of bootstrapping time step $\epsilon$?

We present a table showing how the performance of bootstrapping varies with the time gap $\epsilon = r - t > 0$ between two time steps $r, t$ (on Ising model 5 x 5, $\beta = 0.4$).

As $\epsilon$ decreases, bootstrapping improves performance. However, if $\epsilon$ becomes too small, the scale of the conditional generator $F_{t|r}$ increases rapidly, leading to large fluctuations in the loss magnitude across time. This makes neural network optimization unstable due to inconsistent gradient scales. Hence, selecting a moderately small $\epsilon$ is crucial for stable training.

We appreciate your insightful comment and the connection to existing EBM literature. While these works are relevant, we would like to highlight the key difference in problem setting, as described below.

W1: It is plausible that existing methods from [1, 2, 3] would already be capable of addressing the tasks considered in this paper.

Our task is fundamentally different from the suggested works [1, 2, 3]. The suggested papers focus on energy-based "models", while our work focuses on energy-based "training". The suggested papers aim to train an energy function "from data", while our work aims to train a model "from energy functions".

We hope our responses have addressed your questions and concerns. Please let us know if there is anything we may have missed or if any points require further clarification.

Thank you for following up on our rebuttal and clarifying your point! We do not take it for granted. In what follows, we provide several arguments that show how the suggested energy-based models [1,2,3] do not solve energy-based training (without ground-truth data) and how our contribution is meaningful despite the suggested references.

1. Suggested works [2,3] require initialization from ground-truth data.

The suggested methods [2, 3] critically rely on the assumption of data being available. To be specific, the methods require Markov chains initialized with ground truth samples to sample from the parameterized energy-based models. Without the ground truth samples, which are a good approximation of samples from the energy-based models, the algorithms would require a prohibitive amount of computation. Hence, their performance would degrade significantly for image generation or our settings where ground truth data is unavailable.

Below, we concretize our arguments for each of the suggested methods:

• CDRL [2] trains a generator (the recovery-likelihood model) by running MCMC chains initialized from noise-added ground-truth data.
• DDAEBM [3] trains a variational generator by aligning it to the conditional EBM via a divergence, using noise-added ground-truth data.

2. Cooperative training [1] is vulnerable to mode collapse in our setting.

Cooperative training [1] updates the generator by running short MCMC chains initialized from its own samples and maximizing the MLE objective on the resulting MCMC samples. However, if the generator misses some modes, short MCMC cannot reach them, and the model fails to correct its bias—leading to mode collapse.

We empirically compare cooperative training with our method (EGM + BS) on a continuous 32-dimensional many-well energy function [5]. As shown below, our method discovers more modes and achieves lower energy, indicating better coverage of the energy landscape.

While cooperative training achieves a lower sample $\mathcal{W}_2$ score, this does not reflect comparable performance. The low score results from mode collapse, where samples are concentrated around a single mode—a known failure case for $\mathcal{W}_2$. In contrast, our method achieves broader coverage of the distribution with substantially higher sample quality.

3. Our work is more general. We also remark that our main focus is to develop a general framework that allows sampling both discrete and continuous variables for a broad class of continuous-time Markov chains (diffusion, flow, and jump), while the suggested algorithms are limited to diffusion-based or one-shot generation of continuous variables.

[1] Cooperative training of descriptor and generator networks.

[2] Learning energy-based models by cooperative diffusion recovery likelihood.

[3] Improving adversarial energy-based model via diffusion process.

Thank you for your thoughtful and constructive review. We sincerely appreciate your recognition of the originality of our method. Please find our detailed responses to your comments below.

W1: Energy estimator still has high variance, so why is bootstrapping better?

Thank you for raising this important concern. Our bootstrapping helps training through a trade-off between approximation bias and estimation variance. That is, our algorithm replaces the high variance Monte Carlo approximation of intermediate energy $\mathcal{E}_r$ with a biased surrogate model $\mathcal{E}_r^\phi$. Similar to actor-critic training in reinforcement learning, this trade-off helps stabilize the training of the marginal generator by replacing the high-variance (unbiased) target with a low-variance (biased) target.

To further alleviate your concern, we formally prove that the bootstrapping estimator achieves a lower variance compared to the original estimator.

Proposition. Consider a fully trained surrogate energy model $\mathcal{E}_r^{\phi}(x_r) = \mathbb{E}[\hat{\mathcal{E}}_r(x_r)]$, where $\hat{\mathcal{E}}_r(x_r)$ is a biased energy estimator (from eq. 20). Let $\hat F_t(x_t; \mathcal{E}_r)$ be the bootstrapping estimator using $\mathcal{E}_r$ as the intermediate energy function:

Then, the bias and the variance bound of the bootstrapping estimator with surrogate model $\hat F_t(x_t; \mathcal{E}_r^\phi)$ are given by:

where $\hat F_t(x_t)$ is the estimator without bootstrapping. Since $\text{Var}[\tilde p_r(x_r)] < \text{Var}[\tilde p_1(x_1)]$, the variance of bootstrapping estimator is smaller than the $\hat F_t(x_t)$.

proof sketch) The bias of the energy estimator $\hat{\mathcal{E}}_r$ is given as follows (refer to Corollary 3.2 of BNEM [25]):

Thus, we can approximate the surrogate energy as,

Plugging the equation into the $\hat{F}_t(x_t;\mathcal{E}_r^\phi)$, we get,

Here, $b(x_r^{(i)})$ is close to 0 and concentrated to:

\hat F_t(x_t;\mathcal{E}_r^\phi) \approx (1 - m_b^2)\frac{\sum_i \exp(-\mathcal{E}_r(x_r^{(i)}))F_{t|r}^{x_r^{(i)}}(x_t)}{\sum_i \exp(-\mathcal{E}_r(x_r^{(i)})))} = (1 - m_b^2) \hat F_t(x_t;\mathcal{E}_r)

\text{Bias}[\hat F_t(x_t;\mathcal{E}_r^\phi)] = \mathbb{E}[\hat F_t(x_t;\mathcal{E}_r^\phi)] - F_t(x_t) = \text{Bias}[\hat F_t(x_t;\mathcal{E}_r)] - m_b^2 \left(F_t(x_t) + \text{Bias}[\hat F_t(x_t;\mathcal{E}_r)]\right) = \text{Bias}[\hat F_t(x_t;\mathcal{E}_r)] + O\left(\frac{1}{K^2}\right)

\text{Var}[\hat F_t(x_t;\mathcal{E}_r^\phi)] = (1-m_b^2)^2 \text{Var}[\hat F_t(x_t;\mathcal{E}_r)] \approx \text{Var}[\hat F_t(x_t;\mathcal{E}_r)]

\text{Var}[\hat F_t(x_t;\mathcal{E}_r)]=\frac{4\text{Var}[\tilde p_r(x_r)]}{p_t^2(x_t)K}(1 + \lVert F_t(x_t)\rVert)^2,

\text{Var}[\hat F_t(x_t;\mathcal{E}_r^\phi)] \approx \text{Var}[\hat F_t(x_t;\mathcal{E}_r)] = \frac{\text{Var}[\tilde p_r(x_r)]}{\text{Var}[\tilde p_1(x_1)]}\text{Var}[\hat F_t(x_t)].

Thank you for your reply. On the theory side, unfortunately the reply does not address my main concern.

• Eq. (18) makes the approximation q_{1|r} \propto p_{r|1} for a large jump in time from r to 1 (assuming t is close to r, to allow for small variance in that step).

Previously, in line 135, you stated that approximating q_{1|t} \propto p(t|1) leads to "high variance in the importance weight." So you then make a smaller step to get p_{r|t} which you approximate by p(t|r). When t and r are close, the variance of these new weights are small.

But you are still left with the problem of q_{1|r} \propto p_{r|1}.

You address this problem by introducing bias and reducing variance using (20) - I don't dispute this argument. But it does not address the end-to-end performance of the procedure as a whole, or the reason for introducing the intermediate "r". For example, why not introduce the same eq. 20 bias/variance trick earlier for q_{1|t} \propto p(t|1) ? Why the intermediate step at r?

Re the experiments: thank you for the new experiments - these strengthen the paper.

Due to the continued theory concern, I will maintain my score.

Thank you for following up and restating your concern. Below are detailed responses to your concerns.

Q1. Why introduce the intermediate step $r$ instead of applying the Eq. 20 bias-variance trick directly to generator estimation?

We already apply the analogue of Eq. 20’s bias–variance trade-off to generator estimation via self-normalized importance sampling (SNIS; Eq. 9). However, when temporal gap $t\rightarrow 1$ is large, direct generator estimation remains inaccurate because the estimator depends on two separate sources of variability: the unnormalized density $\tilde p_1(x_1)$ and the conditional generator $F_{t|1}^{x_1}(x_t)$. The combined effect of these two factors leads to high variance that SNIS alone cannot suppress. Bootstrapping addresses this by splitting the high-variance $t \rightarrow 1$ estimation into two lower-variance subproblems, at the cost of a bias from the surrogate energy.

Q2. But you are still left with the problem of $q_{1|r}\propto p_{r|1}$. How does energy estimation with $q_{1|r}$ address the high-variance problem?

As we noted, introducing an intermediate step $r$ through bootstrapping decomposes the high-variance generator estimation from $t\rightarrow 1$ into two subproblems: $t\rightarrow r$ and $r\rightarrow 1$. As you have acknowledged, this effectively reduces the variance in the first step $t\rightarrow r$. For the second step $r\rightarrow 1$ involving energy estimation, we outline several reasons why this estimation problem has lower variance.

1. Energy estimation at time $r$ is easier than the estimation at time $t$.

We set $\epsilon$ to a modest value (e.g., $\epsilon\approx 0.1$) so that $r$ lies closer to $1$ than $t$. This ensures that the $t\rightarrow r$ estimation retains low variance, while the $r\rightarrow 1$ estimation also benefits from a smaller temporal gap. By narrowing this gap relative to $t\rightarrow 1$, the variance of the energy estimator at $r$ is reduced.

2. Energy estimator $\hat{\mathcal{E}}_t$ has lower variance than the marginal generator estimator $\hat F_t$.

Even for the same temporal gap, the energy estimator $\hat{\mathcal{E}}_t$ exhibits lower variance than the generator estimator $\hat F_t$. Energy estimation depends only on unnormalized density $\tilde{p}_1(x_1)$, while generator estimation additionally involves the conditional generator $F\_{t|1}\^{x_1}(x_t)$. This extra source of variability makes the generator estimator more sensitive to the temporal gap. As a result, the variance of $\hat F_t$ is larger than that of $\hat{\mathcal{E}}_t$; we provide a formal argument for this inequality below.

Proposition. Let $\tilde p_1(x_1)$ be an unnormalized density and $F_{t|1}^{x_1}(x_t)$ a conditional generator, evaluated for $x_1 \sim q_{1|t}$. Suppose $\tilde p_1(x_1)$ and $\big\lVert \tilde p_1(x_1)F\_{t|1}\^{x_1}(x_t)\big\rVert$ are sub-Gaussian. Then the variances of the energy estimator and the generator estimator satisfy:

where $K$ is the sample size.

A derivation of the variance for the generator estimator is provided in our response to Reviewer Vwnf (item W3-1), and that for the energy estimator can be found in Corollary 3.2 of [1]. Comparing the leading terms, $\text{Var}[\hat F_t]$ is strictly larger than $\text{Var}[\hat{\mathcal{E}}_t]$ since $4\left(1+\lVert F_t(x_t)\rVert \right)^2>1$.

3. The accuracy of energy estimates can be further improved.

Unlike generator estimation, the accuracy of energy estimation can be enhanced using established techniques. In our method, the forward-looking parametrization evaluates the target energy at intermediate states, yielding an effective approximation to the true intermediate energy. Moreover, as shown in [1], bootstrapping can be applied again to energy estimation, which provably reduces both variance and bias (Proposition 3.4 in [1]).

4. We replace the MC energy estimator with variance to the learned biased surrogate.

To further mitigate variance in energy estimation, we replace the raw Monte Carlo energy estimator with a learned biased surrogate $\mathcal{E}_r^\phi$ within the bootstrapping. As discussed in our previous response, combining bootstrapping with the biased learned surrogate effectively reduces the variance of generator estimation, supporting our claim that the overall bootstrapping framework achieves a lower-variance estimator.

We hope these clarifications address your concerns. If any part remains unclear, please let us know during the remaining discussion period, and we will be glad to elaborate further. Thank you again for your thoughtful feedback.

[1] OuYang et al. BNEM: A boltzmann sampler based on bootstrapped noised energy matching, ICLR 2025.

Thank you for your thoughtful and encouraging review. We appreciate your recognition of our theoretical contributions. Please find our detailed responses below.

W1, Q1: Missing baseline (e.g., iDEM, PIS, DDS, and LEAPS)

We agree on the importance of more baselines; however, we would also like to emphasize that EGM is not directly comparable with the baselines. EGM is a general framework applicable to discrete, continuous, or even discrete-continuous hybrid tasks, while the baselines focus on one of the discrete or continuous tasks.

Nevertheless, we provide a comparison with the baselines. For continuous tasks, we compare with iDEM and PIS on the many-well 32 [6]. We omit DDS from the comparison since it performs similarly to PIS. For discrete tasks, we compare with LEAPS on the 15 x 15 Ising model.

Continuous task (manywell 32):

As the results show, EGM matches or slightly outperforms existing diffusion-based samplers (iDEM, PIS) in the continuous setting.

Discrete task (Ising 15 x 15, $\beta = 0.28$):

For this case, LEAPS is on par with EGM+BS (bootstrapping). However, one can also observe that LEAPS + MCMC is better than EGM+BS at the cost of using more energy calls (45B vs. 85B).

W2-1, Q3: Scalability to high-dimensional problems?

Please refer to our response to W1, Q1, where we present results on a higher-dimensional (d = 225) Ising model to support the scalability of EGM.

Additionally, our primary goal is to expand the design space (e.g., probability path, Markov transition type) and the application of energy-driven training (e.g., multi-modal case), rather than scaling a method to higher dimensions than previous research. Scaling to higher-dimensional tasks remains a challenge for all energy-driven training algorithms.

W2-2, Q2-3: Analysis on training instability is missing.

While the theoretical analyses on stability are outside of the scope of our paper, to alleviate your concern, we provide an empirical analysis of the training stability.

Applying EMA to the energy model can stabilize the training, and we provide a result with and without EMA (on Ising 10 x 10, $\beta=0.4$). All the other hyperparameters are fixed.

Using EMA on the energy-model parameters reduces variance in energy w1 and magnetization score. EMA addresses the moving‑target issue when the generator is regressed on a learned energy model. It is a well‑established technique in reinforcement learning for stabilizing value‑function updates.

W3-1, Q2-2: Lack of bias bounds introduced by self-normalized importance sampling.

We show that the SNIS estimator's error decays as $O(1/\sqrt{K})$, and its bias and variance decay as $O(1/K)$, where $K$ is the sample size. Below, we sketch the proof, and the full rigorous derivation will appear in the revised version.

Proposition. Consider an unnormalized density $\tilde p_1(x_{1})$ and a conditional generator $F_{t|1}^{x_1}(x_{t})$ evaluated on a sample $x_1 \sim q_{1|t}$. Suppose $\tilde p_1(x_1)$ and $\Big\lVert \tilde p_1(x_1) F_{t|1}^{x_1}(x_t) \Big\rVert$ are sub-Gaussian. Then there exists a constant $c(x_t)$ such that with probability at least $1-\delta$,

where $\hat F_t(x_t)$ denotes the SNIS estimator from eq. 9 using $K$ samples $x_1^{(1)}, \\ldots, x_1^{(K)}\sim q_{1|t}$:

Moreover, its bias and variance are expressed as:

proof sketch) By Hoeffding's inequality for sub-Gaussian random variables, there exists a constant $C$ such that

with $1 - \delta$ probability. (Here, $C = \sqrt{\text{Var}[\tilde p_1(x_1)]}$ is a possible choice.)

Since we have bounds for both numerator and denominator of $\hat F_t$, we can also bound the error of $\hat F_t$ to $F_t$. The result is as follows:

Now, let $\hat{A}$ and $\hat{B}$ denote the numerator and the denominator of $\hat F_t$, respectively. We also let $\mathbb{E}[\hat{A}] = A$ and $\mathbb{E}[\hat{B}] = B$. For sufficiently large $K$, the Taylor expansion is expressed as:

To derive the final equation for the bias term, one can express the bias term as follows:

Since $\text{Cov}[\hat{A}, \hat{B}] = \text{Cov}\left[\tilde p_1 F_{t|1}, \tilde p_1\right] / K$ and $\text{Var}[\hat{B}]=\text{Var}[\tilde p_1]/K$, one obtains the conclusion for $\text{Bias}[\hat F_t]$.

To derive the final equation for the variance term, one can combine the sub-Gaussianity of $\hat F_t$ with the error bound on $\left\lVert \hat F_t(x_t) - F_t(x_t) \right\rVert$.

W3-2, Q4: Relationship to iDEM is unclear.

While both iDEM and EGM use importance sampling with energy-based weights, they rely on different score identities with distinct applicability.

iDEM is based on the target score identity (eq. 5 in [37]), which requires access to the target energy gradient and is restricted to the VE probability path. In contrast, EGM adopts the denoising score identity (eq. 2 in [37]), which only requires the conditional generator and applies broadly to general Markov models and state spaces.

The two approaches also differ in variance behavior across time. As shown in [37], combining both estimators can reduce variance in data-driven training. Extending this idea to energy-driven training is a promising direction for future work.

In short, EGM complements and extends iDEM by supporting more general settings.

[37] De Bortoli, Valentin et al. Target Score Matching, arXiv:2402.08667, 2024.

Q2-1: Can you provide a theoretical or empirical analysis of convergence guarantees?

While a formal convergence analysis is beyond the scope of this work, to address your concern, we instead provide an empirical study on the convergence behavior of both EGM and the proposed bootstrapping method.

We report the number of training epochs required for convergence. We define convergence as the point at which the evaluation metric (e.g., energy Wasserstein-1) shows no further significant decrease. (on Ising 10 x 10, $\beta = 0.4$)

Q5: Can you provide a more rigorous analysis of when and why bootstrapping reduces variance?

We prove that bootstrapping reduces variance compared to the non-bootstrapping estimator. We first show that the bias and variance of the bootstrapping estimator with a surrogate energy are close to those of the estimator using the true intermediate energy (i.e., the effect of energy estimation bias is negligible). We then show that bootstrapping with the true energy achieves lower variance than the non-bootstrapping estimator.

This result suggests that, when the surrogate energy is well-trained, bootstrapping effectively reduces the variance of generator estimation.

We state the proposition below. A full proof is provided in our response to Reviewer meoi's comment (W1) and will be included in the revised manuscript.

Proposition. Consider an fully trained surrogate energy model $\mathcal{E}_r^{\phi}(x_r) = \mathbb{E}[\hat{\mathcal{E}}_r(x_r)]$, where $\hat{\mathcal{E}}_r(x_r)$ is a biased energy estimator (from eq. 20). Let $\hat F_t(x_t; \mathcal{E}_r)$ be the bootstrapping estimator using $\mathcal{E}_r$ as the intermediate energy function:

Then, the bias and the variance bound of the bootstrapping estimator with surrogate model $\hat F_t(x_t; \mathcal{E}_r^\phi)$ are given by:

where $\hat F_t(x_t)$ is the estimator without bootstrapping. Since $\text{Var}[\tilde p_r(x_r)] < \text{Var}[\tilde p_1(x_1)]$, the variance of bootstrapping estimator is smaller than the $\hat F_t(x_t)$.

We hope our responses have addressed your questions and concerns. Please let us know if there is anything we may have missed or if any points require further clarification.

We’re grateful that our rebuttal helped clarify your concerns. Thank you for your constructive feedback and positive reassessment.