Energy Matching: Unifying Flow Matching and Energy-Based Models for Generative Modeling

We thank the reviewer for their thoughtful comments and helpful suggestions, which allow us to clarify important aspects of our manuscript. Below, we provide detailed responses addressing the raised weaknesses (W) and questions (Q):

W1: The assumption p_t \approx p_\{eq} (Line 84) may not hold in practice, especially early in the training. The manuscript should elaborate more on the feasibility of such an assumption. Perhaps visualise contributions of the tree components in Eq. 4 depending on $t$. My intuition is that this is why the warmup stage is needed in the first place.

We thank the reviewer for highlighting this important point, as it reveals a missed opportunity in our manuscript to clearly explain the motivation behind the warm-up phase. The statement from line 84 that $p_t$ is close to $p_{eq}$ (in the context of $p_t$ being ‘near data manifold’) indeed may not hold in practice, particularly in the early stages of training. However, line 84 does not imply any guarantees for untrained $V(x)$. Instead, the subsection "First-order optimality conditions" (pp. 3–4) describes properties of the desired solution (such as $p_{data} = p_{eq}$ and Eq. 4) rather than assumptions that must hold throughout training. During training, we explicitly measure how far the current distribution $p_t$ deviates from these optimal conditions and use this imbalance as our learning signal.

We welcome the reviewer's suggestion to visualize how the three forces evolve during training. We will include this visualization in supplementary material.

Early in training, when the potential $V$ is still untrained, contrastive terms provide limited guidance because the negative sample distribution $p_t$ may either be too distant from the data manifold (if initialized in noise) or too close due to poor mixing (if initialized directly on the manifold). This perspective clearly highlights the necessity of the warm-up phase identified by the reviewer. As briefly mentioned around lines 143–144, it is crucial to have negative samples properly distributed around the data manifold to probe low-density regions, thereby constructing energy walls that clearly separate these regions from the high-density data manifold. Once the warm-up phase sufficiently aligns $p_t$ with $p_{data}$, the contrastive terms can more effectively enforce the condition $p_t \approx p_{data} \Leftrightarrow p_t \approx p_{eq}$, actively pulling $p_{eq}$ closer to $p_{data}$.

We thank the reviewer and will explicitly add a sentence at the start of the "First-order optimality conditions" subsection, clarifying that these conditions represent the optimization goal (features of the solution) and highlighting the critical role of the warm-up phase.

W2: The manuscript reports only the FID score for unconditional image generation. Reporting precision and recall [a] would further strengthen the results.

We thank the reviewer for this valuable suggestion. While our results measured by FID significantly surpass existing EBMs (halving the best reported FID scores), we acknowledge that FID alone provides a limited perspective. For proteins, we indeed report additional metrics such as diversity and novelty alongside fidelity measures (see Table 4 in Appendix B, "Posterior sampling"). In the revised manuscript, we will explicitly highlight that, given our FID scores are now on par with state-of-the-art methods, future works on image generation with EBMs should incorporate complementary metrics such as precision and recall [a] for a more complete evaluation.

W3: The proposed training procedure (Algorithm 2) increases training time compared to the standard flow matching due to the sampling required by the contrastive objective. This should be quantified and discussed in ablations (and possibly added to limitations).

While the proposed training procedure indeed adds computational overhead compared to standard flow matching due to the contrastive sampling, we view this overhead as justified given the additional capabilities provided by our approach. As illustrated clearly in Figure 5a, the OT phase alone (without contrastive sampling) achieves results similar to flow matching (same objective, different parameterization), and does not require sampling negative examples. The contrastive objective refines the learned landscape, explicitly enabling functionalities that standard flow-based models cannot achieve directly, such as evaluating log-likelihoods in a single forward pass and performing posterior sampling under arbitrary conditions. Thus, if these functionalities are desired, the additional computation cannot be seen purely as overhead but as necessary to enable these key features.

We fully agree that this computational cost should be clearly quantified and discussed. To address this, we provide the following timing breakdown (per batch on CIFAR-10 using a single NVIDIA R6000):

Additional Timings: EM potential evaluation: 0.148 s; EM velocity generation: 0.317 s

In short, the flow component is roughly 2x slower compared to directly training a velocity field via flow matching due to the additional backward pass. While the CD loss computation itself is approximately 100x slower than the flow loss, CD iterations constitute only about 1–5% of total flow iterations (see Training Details). Using our current negative-sample generation approach (no replay buffers, small fixed step integration), the total training time can increase by up to 12 times in the extreme cases, due to both the curl-free parametrization and negative-sample simulations.

Although this overhead remains negligible for shorter tasks (e.g., protein generation, small images), it could become significant in larger settings. Future improvements could incorporate techniques such as persistent contrastive divergence [E1, E2], replay buffers, or adaptive Langevin steps to significantly reduce generation costs of the negative samples. We omitted these common engineering tricks here in favor of clarity of presentation.

We thank the reviewer for highlighting this point. The overhead described here is characteristic of most traditional EBMs, which typically rely on simulating negative samples during training. Although simulation-based training lost popularity in recent years due to its computational demands and the rise of simulation-free approaches, improvements in hardware capabilities are now bringing simulated negatives back into practical relevance, as they provide complementary signals to simulation-free supervision. Importantly, this overhead is incurred mainly during training and the method remains competitive in terms of inference speed (addressed further in W4). Additionally, Energy Matching can naturally extend to discrete and structured domains, such as molecule generation, where recent work [E3] explores replacing Langevin-based MCMC negative sampling with analytical solutions derived from heat equations on structured spaces.

We will clearly emphasize this context in both the manuscript (Introduction) and appendix (Training Details), making it accessible especially to readers less familiar with EBMs.

• [E1] T. Tieleman and G. Hinton, “Using Fast Weights to Improve Persistent Contrastive Divergence” (ICML 2009)
• [E2] C. Jarzynski, “Nonequilibrium equality for free energy differences”
• [E3] T. Schröder et al., “Energy-Based Modelling for Discrete and Mixed Data via Heat Equations on Structured Spaces” (NeurIPS 2024)

W4: The manuscript parameterises the scalar potential as a deep model. This implies that both forward and backward passes are needed for sampling (as in EBMs). This is significantly slower than in the case of direct parameterisation of vector fields (typicall for flow matching/score matching). The manuscript should discuss this. For example present wallclock time vs FID graph for different baselines.

We fully agree that directly reporting parameter counts and sampling steps serves only as a proxy, and that explicitly comparing wall-clock times against achieved FID scores would offer clearer insights. Indeed, our scalar-potential parameterization requires both forward and backward passes for sampling, similar to traditional EBMs, which roughly doubles (as seen in W3) sampling steps time compared to directly parameterizing vector fields.

We sincerely thank the reviewer for emphasizing this important point and will include a clear wall-clock time vs. FID comparison plot in the revised manuscript. This visualization will clearly compare our method with baseline methods we successfully reproduced. Below, we present the measurements obtained on CIFAR-10 (batch size of 128, NVIDIA R6000), using Euler-Heun (325 steps) and Euler-Maruyama (1000 steps):

• $\dagger$ FID 3.34 achieved by Energy Matching differs from the original submission (FID 3.97) due to adjusted hyperparameters: $\tau^\*=1.0$ (vs $0.9$), $\lambda_{cd}=1\times10^{-3}$ (vs $8\times10^{-3}$), and 325 Langevin steps (vs 300). The architecture is unchanged; updated result in camera-ready.

Thank you for your valuable review and constructive comments, recognizing the strengths of our approach, including its novelty and comprehensive evaluation, and for pointing out areas for improvement regarding the OT solver discussion. Below we address the weaknesses (W) and questions (Q) raised, and we commit to clarifying these points further in the revised version.

W1: The paper fails to state how much the modeling relies on the OT solver. The paper does not provide a clear definition of the OT solver. Is the method employed in the paper the unique OT solver? If not, can authors provide ablations on the OT solver?

Remark 2.1 ("OT solver") on page 4 of our manuscript mentions that we use the "POT" linear-programming (LP) backend (Flamary & Courty, JMLR 2021), which computes an exact 1-to-1 OT plan per minibatch (e.g., size 128). We thank the reviewer for pointing out that this explanation might not be sufficiently clear. This choice is not unique; approximations such as Sinkhorn (with thresholding) or random coupling can also be effective.

Let us (over‑idealise) and model a standardised CIFAR‑10 image as a vector $x \in \mathbb{R}^d$ drawn from $\mathcal{N}(0,I_d)$ with $d = 32 \times 32 \times 3 = 3072$, and pair it with an independent Gaussian‑noise vector $z \sim \mathcal{N}(0,I_d)$. Each coordinate of the difference $(z - x)$ is distributed as $\mathcal{N}(0,2)$, thus the squared Euclidean distance satisfies $\|z - x\|\_2^2 \sim 2\chi\_d^2$, with mean $2d$, standard deviation $\sqrt{8d}$, and relative spread $\frac{\sqrt{8d}}{2d} \approx 0.025$. This thin shell means all entries of the cost matrix $C_{ij} = \|z_i - x_j\|_2^2$ are nearly identical, exemplifying the distance‑concentration effect highlighted by [D0]. Consequently, whether we solve the transport exactly (LP), entropically (Sinkhorn), or even use a random matching, the accumulated optimal‑transport cost should be similar. This has been examined empirically:

• [D1] Table. 5: FID 4.44 (LP) vs 4.46 (random) at 100 steps;
• our new CIFAR-10 runs: 3.97 → 4.00 at 300 steps.
• [D1] Fig. D.2: for batch > 16, solver accuracy saturates on the 2‑moons task.
• [D2] App. C.2: LP and Sinkhorn are indistinguishable unless regularisation is extreme.

We kept POT LP solver for zero hyperparameters, robustness, and negligible cost (see W2). We thank the reviewer and will add a dedicated appendix with the discussion.

• [D0] Aggarwal, et al. “On the Surprising Behavior of Distance Metrics in High Dimensional Space” (ICDT 2001)
• [D1] Tong, et al. “Improving and generalizing flow-based generative models with minibatch optimal transport” (TMLR 2024)
• [D2] Terpin, et al. “Learning diffusion at lightspeed” (NeurIPS 2024)

W2: The paper does not provide an analysis of the computational overhead introduced by the OT solver (e.g., what percentage of the training time per iteration is spent on this step). The OT solver itself can be computationally expensive, especially as the batch size increases. This cost is a critical factor for the scalability of the method, especially when training involves the MCMC sampling.

We acknowledge that the original manuscript did not quantify the computational overhead from the OT solver. Additional experiments show that the OT solver accounts for 1.5% of the training iteration time in Phase 1 (Algorithm 1), dropping to a negligible (0.01%) level in Phase 2 (Algorithm 2), where the computational cost is dominated by negative samples generation.

We thank the reviewer and will include this analysis in the revised manuscript.

W3: The paper forthrightly mentions in its limitations that the method requires extra gradient computation with respect to the input, potentially increasing GPU memory usage. This is an important practical consideration for scalability.

The extra costs buy a key EBM benefit: $V(x)$ supplies explicit curvature and likelihood information, unlocking inverse‑problem solvers and manifold analysis. Flow/score models lack this and instead resort to heavier inference tricks: full‑trajectory differentiation in [D3] or repeated noise/data alternations in  [D4] which come with their own limitations.

• [D3] H. Ben-Hamu, et al., “D-Flow: Differentiating through flows for controlled generation,” (ICLR 2024)
• [D4] M. Mardani, et al., “A variational perspective on solving inverse problems with diffusion models,” (ICLR 2024)

Q1: Figure 6 in Appendix A.2 illustrates that the choice of the temperature-switching parameter (from Eq. 3) can impact FID scores. Could you elaborate on the sensitivity of the model's final performance and training stability to the precise value of and the overall shape/schedule of ? Is the selection of primarily empirical, or are there more principled guidelines for its determination based on dataset characteristics or model behavior?

Thank you for raising this insightful question regarding sensitivity to the temperature-switching parameter ($\tau^*$). The parameter $\tau^*$ marks the activation time ($t>\tau^*$) of the Brownian-motion component in our sampling dynamics; prior to this ($t\le \tau^*$), sampling follows deterministic optimal-transport paths toward high-likelihood regions. As described in Algorithm 3 (Appendix B), conditional sampling terms depend explicitly on the temperature $\varepsilon$, vanishing as $\varepsilon \to 0$. Lowering $\tau^*$ (e.g., 0.8) allows earlier conditioning, potentially aiding conditional tasks, but typically reduces performance in unconditional generation. Conversely, setting $\tau^*$ closer to 1.0 reinforces OT regularization near the data manifold, stabilizing training and enhancing unconditional generation quality. Initially, we presented only a brief ablation of $\tau^*$ (Figure 6, Appendix A.2). Encouraged by your feedback, we've expanded these experiments, individually optimizing hyperparameters (e.g., $\lambda_{cd}$) for different $\tau^*$ values. We now report two hyperparameter sets: Set #1 optimized at $\tau^* =0.9$, and Set #2 optimized specifically at $\tau^* =1.0$.

Our updated CIFAR-10 results are:

Based on these expanded experiments, we improved upon our originally submitted results (FID 3.97 $\to$ 3.34); the updated result (FID 3.34) will be included in the camera-ready version.

Our recommended heuristic is to always start with $\tau^* =1.0$ for unconditional generation tasks, ensuring the best balance of training stability and final performance. If additional conditions are present, one may explore lower values ($\tau^* < 1.0$) later, carefully examining if they provide benefits. We will explicitly incorporate this heuristic rationale into the camera-ready version.

Q2: What is the inference process of this method? Does it need both flow matching part and Langevin MCMC part? If so, would it cost too much time? A standalone, explicit sampling algorithm in the main text would improve the clarity and reproducibility of this paper.

We thank the reviewer for the valuable suggestion; our inference procedure smoothly transitions from deterministic sampling (akin to flow-matching, OT paths and no randomness) to stochastic Langevin MCMC sampling, following a discretized temperature schedule $\epsilon^{(n)}$. Specifically, when the temperature is zero ($\epsilon^{(n)}=0$ for $n \Delta t<\tau^*$), the updates are deterministic (flow-matching), and when the temperature becomes positive ($\epsilon^{(n)}>0$ for $n \Delta t \geq \tau^*$), the updates introduce stochastic noise for modes exploration, becoming Langevin MCMC. This unified scheme is explicitly described in Algorithm 3 (currently in Appendix B.1, which we will move to the main text for clarity). Specifically, at each inference step, the updates are:

• Unconditional update:

  $x_m^{(n+1)} \leftarrow x_m^{(n)} - \nabla_x V_\theta(x_m^{(n)})\Delta t + \eta\sqrt{2\epsilon^{(n)}\Delta t}$

• Conditional update:
  The conditional update is identical in form to the unconditional one but uses a composite potential:

  $V_\theta^{(composit)}(x_m^{(n)}) = V_\theta(x_m^{(n)}) +\frac{\epsilon^{(n)}}{\zeta^2}\|y - A(x_m^{(n)})\|^2 + \frac{\epsilon^{(n)}}{\sigma^2}\sum_{k \neq m} W(x_m^{(n)}, x_k^{(n)})$

Here, briefly:

• $V_\theta(x)$ is the learned scalar potential ($\sim$ negative log-likelihood).
• $A(\cdot)$ denotes a forward measurement operator used in inverse problems.
• $W(\cdot,\cdot)$ represents an optional interaction energy between different samples, e.g. encouraging diversity among generated solutions.
• $\varepsilon^{(n)}$ controls the stochasticity of sampling ($\varepsilon^{(n)}>0$ corresponds to Langevin MCMC).
• $\eta\sim\mathcal{N}(0,I)$ introduces Gaussian noise at each Langevin step.

CIFAR-10 unconditional (batch size of 128, NVIDIA R6000 GPU), using Euler-Heun (325 steps) and Euler-Maruyama (1000 steps):

$^\dagger$ as in Q1: latest improved result for our Energy Matching model (to appear in camera-ready, same architecture).

We thank the reviewer for suggesting this improvement and will explicitly include Algorithm 3 in the main text to enhance clarity and reproducibility.

We thank the reviewer again for the thoughtful feedback and are glad our clarifications were helpful. We will incorporate this discussion into the revised manuscript.

Although Energy Matching is slower than OT-FM due to gradient-based velocity parameterization (requiring both forward and backward passes, an inherent ~2x theoretical slowdown for EBMs), the scaling remains favorable as both passes scale linearly with the number of parameters. Moreover, it remains competitive with diffusion models and significantly outperforms SOTA diffusion-EBM hybrids [C1] in both parameter efficiency (50M vs. 150M) and sample quality (FID 3.34 vs. 3.68).

Energy Matching (and any other EBM method) may not be optimal for raw sampling speed, but its computational cost aligns with recent trends in generative AI towards scaling inference to unlock new capabilities [C2, C3].

If further clarification would be helpful, we are happy to discuss.

[C1] Y. Zhu, et al. “Learning Energy-Based Models by Cooperative Diffusion Recovery Likelihood” (NeurIPS 2023)

[C2] B. Zhang, et al. "Improving Diffusion Inverse Problem Solving with Decoupled Noise Annealing" (CVPR 2025 oral)

[C3] Nanye Ma, et al. "Inference-Time Scaling for Diffusion Models beyond Scaling Denoising Steps" (2025)

We thank the reviewer for their constructive suggestions. We provide point-by-point responses to the raised weaknesses (W) and questions (Q).

W1/Q1: The exposition is very dense and barely understandable for someone outside the OT field. / Please improve the general accessibility of the material to readers outside the OT field.

Thank you for raising this concern. We fully agree that providing intuitive explanations is important, especially for readers less familiar with optimal transport. In particular, we will add a concise paragraph clearly emphasizing the following key insight: a central OT goal in our method is to produce a potential landscape with minimal curvature—effectively a clean, gently sloped "path"—between the noise and data distributions. This smooth slope makes it easier and faster to generate high-quality samples using only a few ODE or SDE integration steps. This capability is essential, as generating accurate, high-quality samples allows us to produce reliable negative examples. These negative examples are then used by the contrastive objective to better capture the geometric structure, or "curvature," of the data manifold close to the actual data points. Thus, away from the data manifold, the OT component naturally regularizes the contrastive objective, enhancing both training stability and final generation quality.

We will incorporate this intuitive perspective clearly in the method section before formally introducing the solution via the JKO scheme. Additionally, we will include a brief summary of essential OT concepts in the appendix for quick reference. We welcome any further suggestions to improve clarity.

W2: The paper title seems overstated. I does not read like a unifying theory of flow matching and energy based models, but rather a combination of the two.

Thank you for this valuable remark. We would like to clarify why we chose this specific title and agree that this deserves a clearer explanation. We will incorporate the following rationale explicitly into the camera-ready version to better communicate why we believe Energy Matching unifies flow matching and energy-based models. Our work is framed entirely in the JKO variational picture: one free‑energy functional is minimised at every step, and the temperature parameter $\varepsilon$ decides which regime we are in. For any positive $\varepsilon$ the updates coincide with the dynamics used to train energy‑based models; as $\varepsilon$ approaches zero the very same updates degenerate to the displacement‑only objective that flow matching employs. We are, to the best of our knowledge, the first to show explicitly that flow matching is the zero‑temperature member of this JKO family, and the title is meant to foreground that link.

Examining Algorithm 3 (p. 15) (also discussed in our response to Q2 from reviewer XTSh) reveals a key limitation of flow matching which operates in the $\varepsilon=0$ limit. Every term that would steer the sampler toward a specific observation is multiplied by $\varepsilon$; setting $\varepsilon$ to zero (precisely the flow‑matching limit) makes those terms disappear. The model is then locked into a single transport path dictated by the prior and cannot produce samples from a conditioned posterior without outside tricks and approximations. Highlighting this drawback, and the way our framework overcomes it, is a novel contribution of the paper.

The optimal‑transport portion of the loss provides positive guidance on where mass should move, while the contrastive-divergence loss contributes a negative signal that discourages excursions into low‑density regions. Rather than merely stitching two concepts together we show that they are part of the same variational picture with complementary losses for different regions of the trajectory.

We agree that explicitly clarifying this rationale will strengthen our contribution and help prevent potential confusion regarding the title. To address your valuable feedback, we will prominently emphasize this unified variational interpretation in the abstract and introduction.

W3: The experiments are limited to smaller scale data (imagenet 32x32) and it's not clear how the methods scales. The authors do acknowledge that in the limitations.

Thank you for highlighting this important point regarding scalability. We agree that our current experiments focus primarily on image datasets of moderate resolution. This scale, however, is a widely studied benchmark, making it well-suited for demonstrating clearly the strengths and relative performance of our approach. Indeed, our method significantly surpasses the previous state-of-the-art FID scores for EBMs by achieving more than a two-fold improvement. We believe this demonstrates the potential to extend our method to higher-dimensional images in future work. Importantly, our approach is also well-suited to problems outside image generation, such as protein design, where typical data dimensionality is substantially lower than even moderate-resolution images. This is a significant and practically relevant domain that benefits directly from the explicit scalar-potential structure inherent in our approach, as demonstrated in our protein-generation experiments. Regarding computational complexity of some of the components:

• Optimal‑transport couplings add only linear $O(d)$ overhead. For very high‑dimensional data there is both theoretical and empirical evidence suggesting dropping the solver and using random couplings that do not cost. (please see response W1 for the reviewer XTSh).

• The full Hessian spectrum with $O(d^3)$ is used exclusively for the Local Intrinsic Dimensionality diagnostic and is neither required for training nor for unconditional/controlled generation; practical deployments can drop this step or approximate it as suggested in the limitation section.

We appreciate the reviewer’s comment and will incorporate a deeper exploration of scalability to larger-scale datasets in our outlook and future research directions.

Q2: The exposition of the results could be clearer. You write that you surpass state-of-the-art EBMs. This is technically true, but somehow omits the fact that other methods in the comparison achieve even better results. I am not in favor of only accepting results that beat state of the art, but you need to argue why your results are still relevant and what advantage/disadvantage other methods have compared to yours.

We thank the reviewer for pointing out that our comparison to non-EBM methods could indeed be clarified further. The inclusion of non-EBM results primarily serves as a contextual benchmark to illustrate that our method significantly advances energy-based modeling toward competitive generation quality, which was previously challenging for EBMs. Within the category of EBMs, our approach clearly achieves state-of-the-art performance. The relevance of EBMs lies in their unique advantages: they inherently provide explicit likelihood estimates (in a single forward pass) for any given sample, enabling straightforward and flexible conditional sampling. Conditions or new interactions between samples can be incorporated by simply adding energies, resulting in a composite energy landscape (explicitly detailed in Algorithm 3 of Appendix B: Posterior Sampling, Section 3.2 on Inverse Problems, and our response to reviewer XTSh's question Q2). Additionally, EBMs naturally provide Hessian-based curvature information about the data manifold, enabling applications (such as LID estimation or more in the response to reviewer 4N5t(Q4)) that flow- or score-based methods cannot directly offer without additional approximations.

However, these advantages come with additional computational overhead during training, which we explicitly acknowledge (see our detailed discussion in W3 of reviewer FNAJ). Whether this overhead is justified depends on the application. For use-cases requiring precise control over generation—such as many scientific applications already outlined in the outlook of our manuscript—this additional cost is typically justified. Similarly, scenarios where conditions, interactions, or labels for conditional generation may not be available at training time, necessitating a model encoding a versatile prior that can accommodate new constraints on demand, further highlight the practical value of our approach.

We sincerely thank the reviewer for emphasizing this important point; we will use the additional page in the camera-ready version to explicitly clarify how to interpret comparisons between EBM and non-EBM methods at the beginning of the results section.

We sincerely thank the reviewer for the feedback and for raising the score. The reviewer's thoughtful comments have helped us improve the manuscript, and we will incorporate the suggested points into the revised version.

Thank you for the thoughtful review and valuable feedback. Below, we provide point-by-point responses addressing the raised weaknesses (W) and questions (Q). We appreciate your comments on the motivation of certain aspects; we agree that this could be clarified further and commit to improving the paper accordingly in the revised version.

W1: Energy matching needs extra overhead compared to flow matching

We agree with the reviewer that Energy Matching introduces additional overhead compared to flow matching, primarily due to the contrastive-divergence training phase (details in Q2/Q3, including complexity and timings). However, this overhead is justified by enabling functionalities beyond unconditional generative tasks. Specifically, the contrastive phase refines the energy landscape to explicitly represent data likelihood, which provides essential advantages including conditional posterior sampling and data manifold curvature information (see detailed motivation in Q2 and Q4). In contrast, flow-based models lose explicit likelihood representation, limiting their applicability in broader probabilistic inference scenarios.

W2: The motivation can be made clearer

We appreciate this feedback. To clarify our motivation, we will highlight explicitly in the introduction that while flow- and diffusion-based models effectively generate high-quality samples, they inherently lack explicit modeling of the unconditional data score. This absence complicates controlling the generation, likelihood estimation, and data manifold curvature analysis (more details in Q4) —limitations that Energy Matching directly addresses. We will further emphasize this motivation clearly in the revision.

Q1: Why is "time-independence" preferred? While it is true that the energy function does not explicitly depend on time, the schedule implicitly plays the role of time.

The data density we ultimately care about is time-independent; thus, constraining the energy to be time-independent eliminates unnecessary degrees of freedom. A single set of weights covers all noise scales, leading to fewer parameters and improved regularization—our model uses 50M parameters compared to 150M in [C1], yet achieves better (lower) FID scores: 3.34 (our latest CIFAR-10 result, to appear in the camera-ready, more in Q3) versus their 3.68. At inference, this approach also removes the need to precisely estimate the noise level, a challenging issue in inverse problems [C2], since estimating time for the schedule $\varepsilon(t)$ is not critical: it is mostly a piecewise constant schedule and inaccuracies in estimating $\varepsilon(t)$ should not significantly impact sample behavior for $t \gg 1$. Thus, parameterizing the potential energy as time-independent provides a minimal representation from a theoretical perspective and robustness in practice.

• [C1] Y. Zhu, et al. “Learning Energy-Based Models by Cooperative Diffusion Recovery Likelihood” (NeurIPS 2023)
• [C2] H. Chung, et al. “Diffusion Posterior Sampling for General Noisy Inverse Problems” (ICLR 2023)

Q2: If the first phase (OT) is accurate, why is the second phase (contrastive objective) even necessary? What is the overhead of the OT solver, as a function of data points? How does it compare to the overhead of distillation?

We thank the reviewer for this important question. Indeed, as shown in Figure 5a, the OT phase alone effectively aligns transport paths, efficiently guiding samples from noise toward the data distribution. However, the resulting potential from this first phase is not directly meaningful in terms of data likelihood—it merely represents a low-curvature landscape that facilitates efficient transport. Crucially, this learned potential after the OT phase does not encode likelihood of the data in any meaningful way. The second phase (CD objective) is therefore necessary to explicitly encode likelihood information into the energy landscape. This is achieved by exploring the data manifold from multiple directions via negative sampling, creating distinct "energy wells" around regions of high data density (Figure 5b). A meaningful potential, proportional to negative log-likelihoods, is the central advantage of energy-based models over flow-based models. This explicit likelihood encoding enables:

• Flexible conditional generation (e.g. inverse problems, interactions), through posterior sampling.
• Explicit (in a single forward pass) log-likelihood estimation.
• Analysis of the data manifold structure (e.g., intrinsic dimensionality estimation).

We agree this motivation could be emphasized more clearly and will explicitly add this discussion to the revised manuscript.

Regarding the overhead of the OT solver as a function of data points and its comparison to distillation: the OT solver we use (POT linear program) has complexity $O(B^2 + d)$, where $B$ is the minibatch size and $d$ is the data dimensionality. In practice, the scaling is mainly driven by the $B$, as $d$ only contributes linearly through the computation of the cost matrix, which is relatively inexpensive. For batches well beyond the size of 128 used here, regularized Sinkhorn methods can reduce complexity to $O(B\log B + d)$. Ablation studies in existing literature [C3, Appendix C.2, Figure 10] indicate minimal differences between these OT methods, provided the regularization does not severely degrade the transport map. The OT solver accounts for 1.5% of the training iteration time in phase 1 (Algorithm 1). This drops to 0.01% during the main training phase (Algorithm 2), after the contrastive-divergence term, which is significantly more expensive due to the simulation-based generation of negative samples, is activated. Training a dedicated distilled model to approximate the global OT map presents a promising direction for achieving sampling with fewer steps. Indeed, precomputing the OT map upfront could eliminate sampling errors that arise when the transport map is supervised only on minibatches. We agree this warrants a deeper discussion and promise to add an additional remark on this topic in the revised version.

Q3: In inference time, the Langevin-type sampling can be slow. How does the performance scale with the number of Langevin steps?

Please refer to Figure 6 in our manuscript, which shows the relationship between FID and the number of Langevin steps (fixed step size $dt = 0.01$). We observe rapid improvement initially, with performance saturating after about 300–325 steps. While Langevin dynamics: $x_m^{(n+1)} \leftarrow x_m^{(n)} - \nabla_x V_\theta(x_m^{(n)}) \Delta t + \eta\sqrt{2\varepsilon^{(n)}\Delta t},\quad\eta\sim\mathcal{N}(0,I)$ can require many steps to reach equilibrium, our training ensures rapid convergence near the data manifold (via the OT flow objective), after which Brownian motion efficiently explores data modes.

Concrete timing results: Unconditional CIFAR-10 sampling (batch size 128, NVIDIA R6000 48GB GPU), using Euler-Heun (325 steps) and Euler-Maruyama (1000 steps):

$^\dagger$ FID 3.34 (lower is better) is our latest CIFAR-10 result using energy matching. The main improvements from the originally submitted result (FID 3.97) are due to adjusted hyperparameters: $\tau$* increased to 1.0, \lambda_{cd} = 1×10^\{-3}, and the Langevin steps increased from 300 to 325.

We thank the reviewer for the important suggestion. We'll highlight our method's competitive sampling efficiency compared to SOTA flow- and score-based baselines in the revision.

Q4: The Hessian of energy, although expensive to compute, contains valuable information. Using Hessian to estimate the LID is a good example. What other things become possible with energy matching that are impossible/challenging with flow matching?

We appreciate this insightful question highlighting a key distinction: energy-based methods explicitly model a scalar field representing data likelihood, unlike flow-based methods. This explicit potential enables probabilistic reasoning tasks (e.g., intrinsic dimensionality estimation, inverse problems, out-of-distribution detection, model selection), which flow-matching inherently lacks. Energy Matching addresses this by learning a potential proportional to negative log-likelihood, supporting likelihood-based computations. Explicit energy modeling enables various Hessian-based probabilistic tasks, e.g.:

• Active-subspace dimension inference [C3]: accelerating inference by identifying data-relevant directions.
• Laplace approximations/Bayes factors [C4]: efficient model selection via marginal likelihood estimation.
• Stochastic-Newton preconditioned MCMC [C5]: accelerated sampling, especially useful in high-dimensional inverse problems.

We thank the reviewer for highlighting this important consideration. While we acknowledge Energy Matching’s higher training cost due to the CD phase (as further elaborated in our response to W3 of reviewer FNAJ), we agree that clearly outlining this trade-off is important. We will emphasize the probabilistic advantages, particularly in likelihood-based reasoning tasks, more explicitly in the revised manuscript.

• [C3] P. Constantine, et al. “Accelerating Markov Chain Monte Carlo with Active Subspaces” (SIAM 2015)
• [C4] E. Daxberger, et al. "Laplace redux-effortless Bayesian deep learning." (NeurIPS 2021)
• [C5] T. Bui-Thanh, et al. “A computational framework for infinite-dimensional Bayesian inverse problems” (SIAM 2015)