Hamiltonian Score Matching and Generative Flows

Authors: peter  holderrieth
yilun  xu
tommi  jaakkola
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:29:39

Abstract: Classical Hamiltonian mechanics has been widely used in machine learning in the form of Hamiltonian Monte Carlo for applications with predetermined force fields. In this paper, we explore the potential of deliberately designing force fields for Hamiltonian systems, introducing Hamiltonian velocity predictors (HVPs) as a core tool for constructing energy-based and generative models. We present two innovations: Hamiltonian Score Matching (HSM), which utilizes score functions to augment data by simulating Hamiltonian trajectories, and Hamiltonian Generative Flows (HGFs), a novel generative model that encompasses diffusion models and OT-flow matching as HGFs with zero force fields. We showcase the extended design space of force fields by introducing Oscillation HGFs, a generative model inspired by harmonic oscillators. Our experiments demonstrate that HSM and HGFs rival leading score-matching and generative modeling techniques. Overall, our work systematically elucidates the synergy between Hamiltonian dynamics, force fields, and generative models, thereby opening new avenues for applications of machine learning in physical sciences and dynamical systems.

Main paper: 
 introduction: 
 Hamiltonian mechanics is a cornerstone of classical physics, providing a powerful framework for analyzing the dynamics of physical systems [33, 4] . The Hamiltonian formalism has been widely applied in machine learning and Bayesian statistics via Hamiltonian Monte Carlo (HMC) [14, 36, 8] . In this setting, the goal is to sample from a probability distribution π whose density π(x) is known up to a normalization factor. In HMC, one interprets ∇ log π(x) as a force function and plugs it into a Hamiltonian ODE to construct a fast-mixing Markov chain exploring the data space quickly. This makes HMC one of the most powerful sampling algorithms to date [36, 23] .In generative modeling, the density π(x) is unknown, only data samples x 1 , . . . , x n ∼ π are given, and the goal is to learn to generate novel samples from π. Current state-of-the-art models are based on diffusion [41, 42, 45, 21] and enjoy widespread success in image generation [39] , molecular generation [11] , and robotics [10] . Diffusion models learn the score function ∇ log π σ (x) for a range of noise scales σ via denoising score matching (DSM) [48] . This enables one to subsequently generate high-quality samples by following a stochastic differential equation [45] .In light of the success of HMC, it is natural to ask whether the Hamiltonian formalism can also improve generative models or provide novel insights into their construction. Previous works have exploited the connection between Hamiltonian physics and generative modeling for specific force fields [13] . However, these works usually consider particular (fixed) force fields and stay within the diffusion framework.More recently, flow-based generative models such as flow matching have enabled scalable training of continuous normalizing flows (CNFs) [31, 32] . These ODE-based models allow to craft first-order ODEs transforming arbitrary distributions from one to another. Inspired by these successes, we consider the Hamiltonian ODE as a Neural ODE [7] . We show that we can marginalize out velocity distributions and then follow backward ODEs that faithfully recover data distributions. Importantly, we provide an associated theorem that holds for any force field. With the growing success of generative models in physical sciences [11, 49, 1] , it is striking that most approaches do not use existing known force fields -sometimes leading to physically implausible results [1] . A framework that allows to reason natively about force fields in the context of generative models would be very promising for such applications [2, 12] . This work aims to build towards such a deeper integration. We explore the intricate relationship between Hamiltonian dynamics, force fields, and generative models. Specifically, we make the following contributions:1. Hamiltonian velocity predictor: We introduce the concept of a Hamiltonian Velocity Predictor (HVP) and show the utility of HVPs for score matching and generative modeling (Section 3). 
hamiltonian score discrepancy (hsd):: 
 We introduce and validate Hamiltonian score discrepancy (HSD), a novel score matching metric based on HVPs and a corresponding score matching method (Section 4).The inspiration of using Hamiltonian dynamics in machine learning comes from considering a data point x ∈ R d as coordinates of an object in R d . Such an object also has a velocity v ∈ R d . Let π be a probability distribution with density function π : R d → R ≥0 . Assuming unit mass, the energy of such an object is given by its Hamiltonian [33, 4] , defined by:EQUATIONwhere U : R d → R is a potential function. The key idea behind using Hamiltonian dynamics in the context of probabilistic modeling and sampling is to set the potential function as negative log-likelihood, i.e., U (x) = -log π(x). One defines the Boltzmann-Gibbs distribution π BG then as:π BG = π ⊗ N (0, I d ), π BG (x, v) = exp(-H(x, v))/Z = π(x)N (v; 0, I d ),i.e. the product distribution of the data distribution and normal distribution. In particular, one can easily draw a sample z from z ∼ π BG by sampling x ∼ π, v ∼ N (0, I d ) and setting z = (x, v).Hamiltonian dynamics describe how an object described by z = (x 0 , v 0 ) evolves over time. It is defined by the ODE [4] EQUATIONEQUATIONi.e. the change of location is the velocity and the change of velocity is the force (here, equals acceleration as we assume unit mass). Let φ : R 2d × R → R 2d , (z, t) → φ t (z) be the corresponding flow, i.e. the function t → φ t (z) is a solution to the above ODE with starting point z.As one would expect from physics, Hamiltonian dynamics φ t preserve the energy of a system, i.e. H(φ t (x, v)) = H(x, v) for all t, x, v (see proof in Appendix B.1). This physical intuition translates into the fact that Hamiltonian dynamics preserve the Boltzmann-Gibbs distribution, i.e.EQUATIONWe include a derivation of this well-known statement in Appendix C.1 as it is so central to this work.The goal of score matching is to learn the score function ∇ log π from data samples x 1 , . . . , x n ∼ π.As the score function naturally appears in the Hamiltonian ODE (see Equation ( 2)), we interpret it as a force function and denote a parameterized score model by F θ : R d → R d . A natural approach to fitting F θ is to minimize the mean squared error between F θ and the true score weighted by their likelihood under π. This leads to the explicit score matching loss [25] given byEQUATION)This loss cannot be minimized directly as one does not have access to ∇ log π and various scorematching methods differ in how they circumvent not having access to ∇ log π (see Appendix A for a detailed overview).A different approach to score matching is to slightly modify the objective by adding Gaussian noise to the data distribution π to get the noisy distribution π σ (x) = N (x; x 0 ; σ 2 I d )π(x 0 )dx 0 [48] . The objective can then be expressed as denoising score matching:EQUATIONfor a constant C σ . Noised data distributions naturally appear in diffusion models, and the denoising score-matching objective, therefore, became the state-of-the-art method to train diffusion models. However, denoising score matching suffers from high variance leading to long training times as well as computationally expensive sampling [44, 52] . In addition, it would be an unreasonable choice for an application where it is important to learn the original data distribution.The evolution is governed by a parameterized force fieldEQUATIONEQUATIONAs we consider Π as part of the definition of a PH-ODE, we write with a slight abuse of notationEQUATIONfor the solution of the ODE assuming z = (x 0 , v 0 ) ∼ Π. We note that while F θ might refer to a θ-parameterized neural network, we can also set it to a fixed known vector field. Both cases will be important.Velocity prediction is all you need. The crucial idea of this work is that one can use PH-ODEs for both score matching and generative modeling by predicting velocities. For this, we use an auxiliary family of functions V ϕ : R d × R → R d , here usually a neural network. Let us consider the following velocity prediction loss:EQUATIONBy minimizing the above loss over ϕ, we train V ϕ to predict the velocity given the location after running the Hamiltonian ODE with starting conditions defined by Π. For a sufficiently rich class of functions V ϕ , the minimizer V ϕ * of the above loss is the expected velocity conditioned on the location (see Appendix D):EQUATION) As we will see, this quantity is "all you need" for generative modeling and score matching. We call V ϕ * the Hamiltonian velocity predictor (HVP). We note that many of flow-based works also implicitly learn to predict velocities [31] , although they do not consider them as separate states. The PH-ODE with F θ preserves the Boltzmann-Gibbs distribution π BG . 
hamiltonian generative flows (hgfs):: 
 We show that the location marginal of a Hamiltonian ODE is generated via the Hamiltonian Velocity Predictor (Section 5). This leads to a novel generative model generalizing diffusion models and flow matching (Section 6).The Hamiltonian ODE has been widely used in Bayesian statistics in the form Hamiltonian Monte Carlo [14, 36] . However, in such settings, it is assumed that one has access to the score ∇ log π (=force), and the goal is to sample from π (by sampling from π BG ). In machine learning tasks, the inverse is true. For such tasks, one has access to data samples x 1 , . . . , x n ∼ π and the goal is to (1) learn ∇ log π (score matching) or (2) create new samples x ∼ π (generative modeling) -or both.Parameterized Hamiltonian ODEs (PH-ODEs). This inspires the definition of a parameterized Hamiltonian ODE (PH-ODE). A PH-ODEs consists of two components:1. Initial distribution Π: The starting condition z = (x, v) ∼ Π is distributed according to a joint location-velocity distribution Π such that its marginal over x is π:EQUATIONIn the default cause, Π equals the Boltzmann-Gibbs distribution π BG , i.e. Π = π ⊗N (0, I d ).The velocity given the location after running the PH-ODE with F θ is zero if starting conditions z = (x 0 , v 0 ) are sampled from π BG :EQUATIONA proof can be found in Appendix C and is based on the fact that test functions linear in v have vanishing expectation if Equation ( 13) holds. We note that the equivalence of conditions (1) and ( 2) is well-known in statistical physics.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes] Justification: The assumptions are fully stated in the supplementary material and referenced. Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
oscillation hgfs:: 
 As special HGFs, we study Oscillation HGFs, a simple generative model rivaling the performance of diffusion models due to in-built scale-invariance (Section 7).2 BackgroundIn this section, we provide a novel score-matching method using PH-ODEs and velocity predictors. We will first connect the score function to a preservation property of Hamiltonian systems (Section 4.1), and then introduce a new score-matching objective derived from this property (Sections 4.2 and Section 4.3).The conservation of the Boltzmann-Gibbs distribution π BG (see Equation ( 4)) is the crucial property that enables Hamiltonian Monte Carlo as it allows for proposals of distant states enabling fast mixing of a Markov chain. The inspiration of this work was to take the inverse perspective: rather than considering the preservation of π BG a useful property of the Hamiltonian ODE, we ask whether it is the defining property of the score function. In other words, is any vector field that preserves the Boltzmann-Gibbs distribution under a PH-ODE automatically the score? And if yes, could we train for this property to learn a score? As we will show, we can characterize the score with an even easier-to-train preservation property solely depending on the velocity predictor. Theorem 1. Let T > 0 and F θ (x) a force field. Let Π = π BG = π ⊗ N (0, I d ). The following statements are equivalent:1. Score vector field: The force field F θ equals the score, i.e. F θ (x) = ∇ x log π(x) for π-almost every x ∈ R d .Condition (3) in Theorem 1 naturally motivates a new way of training F θ to approximate a score. Specifically, our goal is train the force field F θ such that its optimal velocity predictor is zero. By Theorem 1, it necessarily holds F θ = ∇ log π in this case. Unfortunately, such a bilevel optimization is not tractable with stochastic gradient descent in general, as it contains two different objectives.However, a simple trick allows us to convert the above into a single objective. For this, we define the Hamiltonian Score Matching loss:EQUATIONwhereC(θ, t) = E[∥v θ t ∥ 2 ].As the value of C(θ, t) is a constant in ϕ, it holds that the optimal velocity predictor is also the unique minimizer of the Hamiltonian Score Matching loss L hsm (ϕ|θ, t). However, while the argmin is the same, the actual obtained minimum value is drastically different as the next proposition shows. Proposition 1. For a sufficiently rich class of functions (V ϕ ) ϕ∈I , it holds thatEQUATIONThe proof relies on plugging the identity of V ϕ * (see Equation ( 12)) into Equation ( 14)) and can be found in Appendix E. By condition (3) in Theorem 1 (see Equation ( 13)), we want to minimize D hsm (θ|t, π) in order to learn scores. For this, let's define a distribution λ with full support over [0, T ) for T ∈ R >0 ∪ {∞}). With this, we define the Hamiltonian score discrepancy (HSD) asEQUATIONNote that the discrepancy is defined for an arbitrary (regular) vector field F θ -not restricted to scores of probability distributions. By Theorem 1, the discrepancy fulfills all properties that we would expect from a discrepancy to hold: D(θ|π) ≥ 0 for all θ and D(θ|π) = 0 if and only if F θ = ∇ log π. We summarize the findings in the below theorem. The full proof is stated in Appendix F. The Hamiltonian score discrepancy gives a natural measure of how far a vector field is from the desired score vector field. However, at first, it seems rather abstract.The following proposition shows that minimizing this measure has a very intuitive interpretation. In fact, it is closely connected to the explicit score matching loss L esm (see Equation ( 5)). Proposition 2 (Taylor approximation of HSM loss). There exists an error term ϵ(t) such thatEQUATIONand limt→0 1 t 2 |ϵ(t)| → 0.A proof can be found in Appendix G. Intuitively, minimizing the Hamiltonian score discrepancy, therefore, consists of pushing the parabola in Equation ( 18) down onto the x-axis. The above theorem also indicates the optimal choice of T : one should choose T high enough to have a loss value high enough to give signal but low enough to avoid errors due to numerical integration of the ODE.Beyond its theoretical value, we can explicitly minimize the HSD, a method we coin Hamiltonian Score Matching (HSM). To minimize the HSD, two networks V ϕ and F θ can jointly optimize Equation (14) . There are two difficulties coming along with this: (1) One has to simulate trajectories. This can be done via Neural ODEs [7] with constant memory. (2) One has to run a min-max optimization. Here, a big toolbox developed for GANs for training stabilization can be used [35, 19] .On the other hand, we hypothesize that HSM has two advantages: (1) every trajectory of HSM gives several points of supervision effectively augmenting our data and (2) we can learn the original ("unnoised") data distribution π. However, please note that we do not propose Hamiltonian Score Matching as a replacement for denoising score matching in diffusion models. Rather, it is a scalable alternative to score matching methods that learn the original ("unnoised") data distribution π.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: We provide experimental details in the supplementary material (Appendix L). Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.• Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [No] Justification: The core contribution of our paper is a method. The extent of the societal impact via deepfakes is not changed due to this work. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
hamiltonian generative flows: 
 Next, we show that training a general velocity predictor of a Hamiltonian ODE is useful even if F θ ̸ = ∇ log π. This leads to a generative model that we coin Hamiltonian Generative Flows (HGFs). As F θ = F is fixed and not trained here, we write (x t , v t ) = (x θ t , v θ t ) for the solution of the PH-ODE. Let us denote Π(x, v, t) as the distribution of (x t , v t ) at time t and the location marginalEQUATIONThe location marginal describes a probability flow starting from our data distribution π = π(•, 0). It turns out that the optimal velocity predictor is exactly the vector field that generates π(x, t).Proposition 3. Let π T be the distribution such that x θ T ∼ π T . Let V * ϕ be the Hamiltonian Velocity Predictor (see Equation (12) ). Then by sampling x T ∼ π T and running the velocity predictor ODEEQUATION)backwards in time, we will have x(0) ∼ π, i.e. we can sample from the data distribution π. More specifically, the optimal velocity predictor V ϕ * generates the probability path π(•, t).The proof uses the fact that the vector field G(x, v) = (v, F θ (x, t)) is divergence-free to show that V ϕ * fulfils the deterministic Fokker-Planck equation (see Appendix H). The above proposition allows us to build a generative model by training an HVP. We coin this model Hamiltonian Generative Flows (HGFs). To make this framework tractable, we need two criteria to be fulfilled:1. (C1) Forward ODE efficiently computed: For efficient training, we need to be able to compute x θ t , v θ t efficiently -either via an analytical expression or ODE solvers. 6 Diffusion Models and Flow Matching as HGFs with zero force field 
diffusion models as hgfs: 
 We can recover diffusion models with a variance-preserving (VP-SDE) noising process [45] as a special case of HGFs. If we simply set Π = π BG and F θ (x) = 0 -no force applied. In this case, we get that x t = x + tv and v t = v leading to the training objective:E x∼π,v∼N (0,I d ) [∥V ϕ (x t , t) -v t ∥ 2 ] =E x∼π,ϵ∼N (0,I d ) [∥V ϕ (x + tϵ, t) -ϵ∥ 2 ]which equals denoising score matching (Equation ( 6)) with score model -tV ϕ (x, t) = ∇ log π σ(t) with σ(t) = t. In this case, Hamiltonian Generative Flows correspond to training a diffusion model and the velocity predictor corresponds to a denoising network (often denoted as ϵ θ in DDPMs [21] ). It then holds x T ∼ π T ≈ N (0, σ 2 (t)I d ) and the velocity predictor ODE then reduces to the well-known probability flow ODE formulation of diffusion models with noise schedule σ(t) = t:x T ∼ π T ≈ N (0, σ 2 (t)I d ), d dt x(t) = -σ(t)σ(t)∇ x log π σ(t) (x) = V (x(t), t)In fact, the above is a universal way of modeling diffusion models [26] . In other words, diffusion models are a special case of HGFs for the zero-force field. In this perspective, different diffusion models correspond to different time rescaling and preconditioning of the network. The location marginals π(x, t) only fully converge to a Gaussian in the limit of t → ∞ (see Figure 2 ).Flow matching with the CondOT probability path [31] , a current state-of-the-art generative model, can be easily considered an HGF model. As in diffusion, let us consider the zero force field F θ and let's consider a coupled initial distribution Π:EQUATIONSimilarly, the velocity prediction loss corresponds to the OT-flow matching loss:EQUATIONand flow model corresponds to the velocity predictor ODE. Therefore, diffusion models and OT-flow matching are both HGFs with the zero force field -the difference lies in a coupled construction of the initial distribution (see Figure 2 ). The coupled construction allows OT-flow matching to have straighter probability paths, leading to improved generation quality for the same number of steps [31] . 
oscillation hgfs: 
 So far, we studied optimal velocity predictors V ϕ for two extreme cases: either F θ = 0 or F θ = ∇ log π. Finally, we want to investigate a different choice of F θ to construct HGFs. Here, we study Oscillation HGFs that correspond to a natural extension. In Appendix I, we give another example that we coin Reflection HGFs.A simple design of a force field is to use the physical model of a harmonic oscillator, i.e. to set F θ (x) = -α 2 x with Π = π BG and α > 0. The flow of the ODE then becomes:EQUATIONI.e. condition (C1) is fulfilled as we can simply compute the ODE analytically. Setting T = π/(2α), it holds that (x t , v t ) = (v, -αx). In particular, π T = N (0, I d /α 2 ) -condition (C1) is easily fulfilled. Therefore, the above choice gives us a natural generative model based on harmonic oscillators that we coin Oscillation HGFs. To summarize, they have the following simple training objective:EQUATIONA natural choice for α is to setα = d/E x∼π [∥x∥ 2 ]. With this, the scale of the n-th derivative (including n = 0) of the inputs and outputs in the training objective is constant in time (see Figure 2 ), i.e. for all t = 0:EQUATIONIn the context of critically-damped Langevin diffusion [13] , it was already observed that a constant scale in velocity space leads to improved training and better performance. Here, we extend this idea of a constant scale from the velocity to the n-th derivative. 
related work: 
 Assessing and training energy-based models. Stein's discrepancy [16] is a well-known measure to assess the quality of energy-based models based on Stein's identity [46] . Based on this metric, [17] developed a method that is similar to ours where a critic is optimized to assess the quality of an energy-based model via Stein's discrepancy and jointly trained with the energy model via min-max optimization. [25] introduced score matching as a method by showing that the explicit score matching loss (see Equation ( 5)) can be implicitly trained if one computes the trace of Hessian of the energy function -an expensive step. To expedite this, [34] introduced curvature propagation for an unbiased Hessian estimate, while [43] used Hutchinson's Trick to estimate the trace. In practice, both methods suffer from high variance due to their underlying Monte Carlo estimators.Flow matching and Stochastic interpolants. As already seen for a special case in Section 6, HGFs are strongly connected to Flow Matching [31] and Stochastic Interpolants [3] . They construct probability paths that fulfill the continuity equation by predicting derivatives of flows (i.e. velocities) in the same way how in this work, we predict velocities as marginals of an extended state space. The differences of these 3 works lie in the design perspective: HGFs consider 2nd-order ODEs in an extended state space R d × R d with a simple initial velocity distribution (here, N (0, I d )), while FM considers 1st-order ODE paths in R d converging to a simple final location distribution. Flow matching conditions on final states (usually at t = 1), while our framework conditions on the velocity of the current state (arbitrary t) and is centered around forces. This work arrives at the ideas of conditional velocity predictors via the search of properties that are conserved under Hamiltonian dynamics (see Theorem 1). We elucidate the mathematical connection in more detail in Appendix J.Generative models and Hamiltonian physics. V-Flows also consider augmenting the state space with velocities deriving an ELBO objective for CNFs [6] . [13] extend diffusion models to joint state-velocity samples that converge to a joint normal distribution. One difference is that we only need to run the backward equation in state space R d as opposed to extended state-velocity space R d × R d . Though rather unmotivated, Oscillation HGFs could, in principle, also be derived as an EDM model with preconditioning [26] . Finally, several works have, like this, explored generative models based on specific physical processes, e.g. Poisson flow generative models [50, 51] . A few works also combined Hamiltonian physics with deep learning. For example, [18] use conservation of energy as an implicit bias to learn networks for physical data. Conversely, deep learning was also used to accelerate HMC sampling [15] , e.g. by training custom MCMC kernels [30] or correct for complex geometries via flows [22] . Very recently, score matching approaches were also designed to leverage existing force fields as part of a diffusion model that samples from an energy landscape [2, 12] .Acceleration Generative Model (AGM). The AGM model [9] also uses constructions in phase space (joint position and velocity space) and 2nd order ODEs. While AGM focuses on learning the force field, our approach primarily focuses on learning the optimal velocity predictor. While we also consider optimizing the force field by minimizing the norm of the optimal velocity predictor, this happens in the "outer loop" of the maximization -the inner loop optimizes the optimal velocity predictor. Further, ATM focuses on bridging two desired distributions by posing a stochastic bridge problem in phase space. We do not consider the problem of bridging distributions. In contrast, our framework centers around energy preservation and divergence from that preservation (for optimal velocity predictors that are not zero). Specifically, we establish a connection to Hamiltonian physics and a property of the preservation of energy. This allows us to introduce a further bi-level optimization and the possibility of joint training for score matching. 
hamiltonian score discrepancy: 
 As we introduced Hamiltonian score discrepancy as a novel score-matching metric, we first empirically investigate our theoretical insights on Gaussian mixtures (see Figure 1 ). As one can see in Figure 3a , the Hamiltonian score discrepancy is highly correlated with the explicit score matching loss. Further, we can validate empirically that the Taylor approximation derived in Proposition 2 is pretty accurate for large t (see Appendix I). Overall, these results indicate that the Hamiltonian score discrepancy is a natural metric to assess score approximations.Further, we investigate whether explicitly minimizing the Hamiltonian score discrepancy leads to accurate score approximations. We jointly train velocity predictors and score networks as described in Section 4. As one can see visually in Figure 1 , this approach can faithfully learn score vector fields. In addition, we investigate the signal-to-noise ratio for gradient estimation. As shown in Figure 3c , the gradient estimates of HSM have significantly lower variance compared to denoising score matching at lower noise levels σ. The reason for that is that we allow for supervision across a full trajectory at locations for the same data points -effectively acting as data augmentation. Of course, this comes at the expense of simulating the Hamiltonian trajectories for ∼ 5 steps. In the form of diffusion models and flow matching, HGFs have already been extensively optimized and achieved state-of-the-art results. Instead, we investigate whether also other non-zero force fields, specifically Oscillation HGFs, can lead to generative models of high quality.We focus on image generation benchmarks. Specifically, we train a Oscillation HGF on CIFAR-10 unconditional and conditional. As two central benchmarks, we use the original SDE formulation of diffusion models [45] as well as the EDM framework [26] , a highly-tuned optimization of diffusion models. Our hypothesis is that Oscillation HGFs should work well out-of-the-box, as the scale of their inputs and outputs stay around constant (c.f. Eq. 25). Therefore, we remove any preconditioning optimized for diffusion models (scaling of inputs and outputs and skip connections) [26] and train on the default DDPM architecture [21] (see details in Appendix L). Our results are encouraging: without hyperparameter tuning, Oscillations HGFs can sample high-quality images and surpass most previous methods (see Table 1 ) measured by Frechet Inception Distance (FID) [20] . While they still lack behind the EDM model, this difference might well be explained by the fact that architectures and hyperparameters have been optimized for diffusion models over several works that are hard to replicate.To investigate whether similar results can be achieved similar performance at higher resolutions, we perform another benchmark on the FFHQ dataset at 64x64 resolution. Here, our results are similar: Oscillation HGFs improve upon original diffusion models with a small performance margin to the EDM model. They can generate high-quality faces that appear realistic (see Figure 4 ). 
conclusion: 
 Our work systematically elucidates the synergy between Hamiltonian dynamics, force fields, and generative models -extending and giving a new perspective on many known generative models. We believe that this opens up new avenues for applications of machine learning in physical sciences and dynamical systems. However, several limitations remain. Minimizing the Hamiltonian Score Discrepancy (HSD) via a default min-max algorithm is scalable but requires adversarial optimization. Future work can focus on adapting the HSD framework, e.g. to develop denoising Hamiltonian score matching that could allow for guaranteed convergence. For HGFs, the extended design space presented has only been explored for data without known force fields (image generation, here). Future work can focus on specific applications that require domain-specific force fields, e.g. for molecular data. Further adaptions might be required in such settings as such data often lie on manifolds. A further challenge is that HGFs not necessarily converge to a known distribution for more complex force fields. Therefore, we anticipate that future work will focus on adapting HGFs and related models to this challenge to design domain-specific models. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [No] Justification: We do not release any data. We will release models upon publication. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? 
 
 Appendix: To make score matching tractable, one can express the explicit score matching loss via [25] LEQUATIONwhere ∇ • F θ is the divergence of the vector field and the constant C is independent of θ. While this loss has been widely used [29] , this objective is hard to optimize with neural networks as the divergence requires to backpropagate d times. Still, the divergence can be approximated via Hutchinson's trick [24] leading to sliced score-matching [43] .In this section, we prove the fundamental properties of the flow φ t of the Hamiltonian ODE (see Equation ( 2)). As these properties are used throughout this work and usually presented in the context of the physics literature, we include the proofs here for completeness, solely expressing it in the language of probability theory. Throughout this section, let J ∈ R 2d×2d be the matrix defined by:J = 0 1 d -1 d 0 B.1 Preservation of energy. Proposition 4. For all t ∈ R it holds that H • φ t = H Proof. We follow [5, theorem 2.2]. For all z ∈ R 2d one has ⟨z, Jz⟩ = 0. Hence d dt H(φ t (z)) = ⟨∇H(φ t (z)), d dt φ t (z)⟩ = ⟨∇H(φ t (z)), J∇H(φ t (z))⟩ = 0which implies that H(φ t (z)) = H(φ 0 (z)) = H(z).Proposition 5 (Symplectic property). For all z = (x, v) ∈ R 2d and t ∈ R the Jacobian Dφ t (z) satisfies the following equation:EQUATIONIn particular, |detDφ t | ≡ 1, i.e. the φ t is volume-preserving.Proof. Here, we follow [33, proposition 2.6.2]. First of all, it can be easily seen that Equation ( 27) is equivalent to the statement that the bilinear form β(u, v) = ⟨u, J T v⟩ on R 2d fulfilsEQUATIONBy taking derivatives of the right-hand side and afterwards using the identities β(Jx, y) = ⟨J 2 x, y⟩ = -⟨x, y⟩ and β(x, Jy) = ⟨x, J T Jy⟩ = ⟨x, y⟩, one gets:d dt β(Dφ t (z)u, Dφ t (z)v) =β( d dt Dφ t (z)u, Dφ t (z)v) + β(Dφ t (z)u, d dt Dφ t (z)v) =β(D d dt φ t (z)u, Dφ t (z)v) + β(Dφ t (z)u, D d dt φ t (z)v) =β(J∇ 2 H(φ t (z))Dφ t (z)u, Dφ t (z)v) + β(Dφ t (z)u, J∇ 2 H(φ t (z))Dφ t (z)v) = -⟨∇ 2 H(φ t (z))Dφ t (z)u, Dφ t (z)v⟩ + ⟨Dφ t (z)u, ∇ 2 H(φ t (z))Dφ t (z)v⟩ =0where in the last step I have used that the Hessian is symmetric. One can concludeβ(Dφ t (z)u, Dφ t (z)v) = β(Dφ 0 (z)u, Dφ 0 (z)v) = β(u, v)Since |det J| = 1, one immediately gets |detDφ t | ≡ 1.This section gives proof that the Hamiltonian ODE defined with the score preserves the Boltzmann-Gibbs distribution. This implication is well-known, and the proof is included for completeness following [5] .For all t ∈ R and Borel sets D, it holds by Propositions 4 and 5 and a change of variables:EQUATIONEQUATIONEQUATIONNote that the symplectic property is crucial, e.g. consider the simple pendulum (U (x) = 1 2 x 2 ) and g(z) = (∥z∥, 0). This function fulfills π BG (g(z)) = π BG (z) but it does not leave the distribution invariant (it is not symplectic, either).We write φ t (z) = (x ∇U t , v ∇U t ) for the solution with force field -∇U . Let's pick an arbitrary sufficiently regular test function f : R 2d → R and let's defineF θ (t) = E z∼π BG [f (φ θ t (z))] F ∇U (t) = E z∼π BG [f (φ t (z))]As the Hamiltonian dynamics with force network -∇U preserve the Boltzmann-Gibbs distribution, F ∇U must be a constant function, i.e. it derivative is zero. Therefore, we can compute:0 = d dt F ∇U (t) = d dt E z∼π BG [f (φ t (z))] =E z∼π BG d dt f (φ t (z)) =E z∼π BG ∇ z f (φ t (z) T d dt φ t (z) =E z∼π BG ∇ z f (φ t (z)) T v ∇U t (z) -∇U (x ∇U t (z))I.e. at time t = 0:EQUATIONAs φ θ t also preserves the Boltzmann-Gibbs distribution, we can follow the same computation to also get that:EQUATIONSubstracting Equation (33) from Equation (34), we get that:EQUATIONwhere DF θ denotes the Jacobian and ∇ 2 U the Hessian. Then inserting Equation ( 36) into Equation (35) we get that:0 =E z=(x,v)∼π BG v T (DF θ (x) + ∇ 2 U (x)) F θ (x) + ∇U (x) T 0 F θ (x) + ∇U (x) =E z=(x,v)∼π BG ∥F θ (x) + ∇U (x)∥ 2 =E x∼π ∥F θ (x) -[-∇U (x)]∥ 2This implies that F θ (x) = -∇U (x) = ∇ log π for π-almost every x.Finally, we show that condition (3) is equivalent to condition (1) . Note that if (x θ t , v θ t ) ∼ π BG = π ⊗ N (0, I d ), then v θ t is independent of x θ t and it holds that E[v θ t |x θ t ] = 0. Therefore, condition (2) trivially implies condition (3). As we have already have shown that condition ( 2) is equivalent to condition (1), it is sufficient to prove that (3) implies (1) . The proof is similar to the proof for (2)⇒(1).Let's assume that the following condition holds for z ∼ π BGEQUATIONwhere for brevity we write(x θ t , v θ t ) = φ θ t (z). Then for f (z) = f (x, v) = v T (F θ (x) + ∇U (x))as selected above, we can equally derive:F θ (t) =E z∼π BG [f (φ θ t (z))] =E[f (x θ t , v θ t )] =E[(v θ t ) T (F θ (x θ t ) + ∇U (x θ t ))] =E[E[(v θ t ) T (F θ (x) + ∇U (x))|x = x θ t ]] =E[E[v θ t |x θ t ] T (F θ (x θ t ) + ∇U (x θ t ))] =E 0 T (F θ (x θ t ) -∇U (x θ t )) =E[0] =0In the same way as above, we can now deduce that d dt F θ (t) |t=0 = 0. Similarly to above, we can complete the proof by using Equation (36) :0 = d dt E z∼π BG [f (φ θ t (z))] - d dt E z∼π BG [f (φ t (z))] =E z=(x,v)∼π BG ∇ z f (z) T 0 F θ (x) + ∇U (x) =E x∼π ∥F θ (x) -[-∇U (x)]∥ 2Again, we can deduce that F θ (x) = -∇U (x) = ∇ log π(x) for π-almost every x.We formally proof the identity in Equation ( 12). We will use the following well-known characterization of the expectation: Lemma 1. Let X ∈ R d be a random variable. Then:EQUATION)Proof. Let c = E[X] ∈ R d , then E[∥X -m∥ 2 ] =E[∥X -c + c -m∥ 2 ] =E[∥X -c∥ 2 + 2(X -c) T (X -m) + ∥c -m∥ 2 ] =E[∥X -c∥ 2 ] + 2(E[X] -c) T (X -m)] + ∥c -m∥ 2 =E[∥X -c∥ 2 ] + ∥c -m∥ 2 ≥E[∥X -c∥ 2 ]This implies the statement.Therefore, we can apply Lemma 1 conditionally on x θ t to see thatV ϕ * (x, t) = E[v θ t |x θ t = x]We can derive:EQUATION)EQUATIONAsE z∼π BG ∥v θ t ∥ 2 is a constant in ϕ, we can see that arg min ϕ L hsm (ϕ|θ, t) = arg min ϕ E z∼π BG [∥V ϕ (x θ t , t) -v θ t ∥ 2 ]Hence, the Hamiltonian velocity predictor is the minimizer of the above objective. Inserting Equation (12) into Equation ( 39), we getmax ϕ L hsm (ϕ|θ, t) =L hsm (ϕ * |θ, t) =E z∼π BG 2E[v θ t |x θ t ] T v θ t -∥E[v θ t |x θ t ]∥ 2 =E z∼π BG 2E[E[v θ t |x θ t ] T v θ t |x θ t ] -∥E[v θ t |x θ t ]∥ 2 =E z∼π BG 2E[v θ t |x θ t ] T E[v θ t |x θ t ] -∥E[v θ t |x θ t ]∥ 2 =E z∼π BG 2∥E[v θ t |x θ t ]∥ 2 -∥E[v θ t |x θ t ]∥ 2 =E z∼π BG ∥E[v θ t |x θ t ]∥ 2where have used that the conditional expectation E[v θ t |x θ t ] is a constant conditioned on x θ t . This finishes the proof.As λ is a distribution with full support over [0, T ), it holds thatD hsm (θ|π) = E t∼λ,z∼π BG ∥E[v θ t |x θ t ]∥ 2 = 0 if and only if for very 0 ≤ t < T (up to measure zero) E[v θ t |x θ t ] = 0 By Theorem 1, this is equivalent to F θ = ∇ log π. Hence, D hsm (θ|π) = 0 if and only if F θ = ∇ log π.Remark. Technically speaking, we maximized V ϕ for a fixed t in Proposition 1. However, we remark that maximizing it across t leads to the same result for all 0 ≤ t < T under reasonable regularity conditions. More specifically, assuming that π is a smooth density with full support, the map (x, t) → E[v θ t |x θ t = x] is continuous in t and x. Therefore, as long as (V ϕ ) ϕ∈I covers all continuous function in t and x, the result above is the same.The proof relies on a Taylor approximation of L hsm (θ, t) around t = 0. To finish the proof, we have to show the following three equations:EQUATION)EQUATION)Proof of Equation ( 43) Note that at time t = 0, it holds that v θ t = v ∼ N (0, I d ) andx θ t = x ∼ π. Therefore, E[v θ t |x θ t ] |t=0 = E x∼π,v∼N (0,I d [v|x] = E v∼N (0,I d ) [v] = 0 Therefore, by Proposition 1 D hsm (θ|0, θ) = E[∥E[v θ t |x θ t ] |t=0 ∥ 2 ] = E[∥0∥ 2 ] = 0Proof of Equation ( 44) We can compute:d dt D hsm (θ|t, π) =E[ d dt ∥E[v θ t |x θ t ]∥ 2 ] = 2E[E[v θ t |x θ t ] T d dt E[v θ t |x θ t ]]So at time t = 0:d dt D hsm (θ|t, π) |t=0 =2E[E[v|x] T d dt E[v θ t |x θ t ] |t=0 ] = 2E[0 T d dt E[v θ t |x θ t ] |t=0 ] = 0Proof of Equation ( 45) Let's take the second derivative:d 2 d 2 t D hsm (θ|t, π) =2 d dt E[E[v θ t |x θ t ] T d dt E[v θ t |x θ t ]] =2E[E[v θ t |x θ t ] T d 2 d 2 t E[v θ t |x θ t ]] + 2E[∥ d dt E[v θ t |x θ t ]∥ 2 ]by the product rule. And at time t = 0:EQUATIONLet π θ t : R d × R d → R be the density of (x θ t , v θ t ) and compute:EQUATION)EQUATION)Let G(x, v) = (v, F θ (x)) T be the Hamiltonian vector field. By the deterministic Fokker-Planck equation [37] , we can derive that:d dt π θ t = -∇ • [π θ t G] = -G T ∇π θ t -π θ t ∇GNote that the Jacobian of F is given byDG(x, v) = 0 I d DF θ (x) 0In particular, G is divergence-free, i.e. ∇ • G = tr(DG) = 0. Therefore,d dt π θ t = -G T ∇π θ t and for particular x, v ∈ R d : d dt π θ t (x, v) |t=0 = - v F θ (x) T ∇π θ 0 (50)As π θ 0 = π BG by construction, we can derive that∇π θ 0 (x, v) =∇π BG (x, v) = 1 Z ∇[exp(-U (x) - 1 2 ∥v∥ 2 )] = - 1 Z ∇U (x) v exp(-U (x) - 1 2 ∥v∥ 2 ) = - ∇U (x) v π BG (x, v)Inserting this into Equation ( 50), we getd dt π θ t (x, v) = v F θ (x) T ∇U (x) v π BG (x, v) =v T (F θ (x) + ∇U (x))π BG (x, v)And hence,d dt π θ t (x) = d dt π θ t (x, v)dv = v T (F θ (x) + ∇U (x))π BG (x, v)dv = (F θ (x) + ∇U (x))π(x)dv T vN (v; 0, I d )dv =0We can insert these identities into Equation (49) to get:d dt E[v θ t |x θ t = x] |t=0 = v π θ t (x) d dt π θ t (x, v) -π θ t (x, v) d dt π θ t (x) (π θ t (x)) 2 dv |t=0 = v π(x)v T (F θ (x) + ∇U (x))π BG (x, v) -0 π(x) 2 dv = v v T (F θ (x) + ∇U (x))π(x) 2 π(x) 2 N (v; 0I d )dv = vv T N (v; 0I d )dv (F θ (x) + ∇U (x)) =I d (F θ (x) + ∇U (x)) =F θ (x) + ∇U (x)Combining this with Equation ( 46), we get that:EQUATIONTaylor approximation Finally, we can combine the above derivations to get a Taylor approximation of L hsm (θ, t) around t = 0, i.e. for ϵ : R → R with limt→0 1 t 2 |ϵ(t)| = 0 we get D hsm (θ|t, θ) = D hsm (θ|0, θ) + t d dt D hsm (θ|t, θ) |t=0 + 1 2 t 2 D hsm (θ|t, θ) |t=0 + ϵ(t) = 1 2 t 2 4L esm (θ; π) + ϵ(t) = 2t 2 L esm (θ; π) + ϵ(t)This finishes the proof.Again, let's consider a probability distribution π : R d → R and the ODE( d dt x(t), d dt v(t)) T = (v(t), F θ (x(t), t))where we know allow F θ to be time-dependent. Let (x t , v t ) be a solution to the above with ODE with(x 0 , v 0 ) = (x, v) ∼ π ⊗ N (0, I d ).In addition, write Π(x, v, t) for the distribution at time t (i.e.(x t , v t ) ∼ Π(•, •, t))and the location marginalΠ(x, v, t)dv = π(x, t)Finally, we writeV (x, t) = E[v t |x t ] = vπ(v|x, t)dv = V ϕ * (x, t)for the optimal velocity predictor.Deriving marginal ODE. We now show that the evolution of the first marginal can be replicated by an ODE that only depends on V . By the Fokker-Planck equation, we can derive forG(x, v) = (v, F θ (x)) T : ∂ ∂t Π(x, v, t) = -∇ x,v • [ΠG](x, v, t) = -G(x, v, t) T ∇ x,v Π(x, v, t) -[∇ x,v • G(x, v, t)]Π(x, v, t) = - v F θ (x) T ∇ x,v Π(x, v, t) = -[v T ∇ x Π(x, v, t) + F θ (x) T ∇ v Π(x, v, t)]where we used in the third equation that ∇ x,v • G = 0 (i.e. that G is divergence-free). Therefore, we can derive that:∂ ∂t π(x, t) = ∂ ∂t Π(x, v, t)dv = -[v T ∇ x Π(x, v, t) + F θ (x) T ∇ v Π(x, v, t)]dv = -v T ∇ x Π(x, v, t)dv -F θ (x) T ∇ v Π(x, v, t)dv = -v T ∇ x Π(x, v, t)dv -F θ (x) T 0 = -v T ∇ x Π(x, v, t)dv = -∇ x • vΠ(x, v, t)dv = -∇ x • [Π(x, t) v Π(x, v, t) Π(x, t) dv] = -∇ x • [Π(x, t) vΠ(v|x, t)dv] = -∇ x • [Π(x, t)V (x, t)]In other words, the vector field V (x, t) satisfies the continuity equation. Therefore, if we initialize x T ∼ π(•, T ) and evolve the ODEd dt x(t) = V (x(t), t)backwards until t = 0, we know that x(0) ∼ π(•, 0) = π -we sample from our data distribution.We briefly show that HGFs can give us models other than Oscillation HGFs. For this, we introduce Reflection HGFs here and plot the result of training them in Figure 5 and Figure 6 . The idea of the model is that particles can move freely in a box without collision with walls ("very strong forces") at the boundaries of the data domain making the particles bounce back (this can be made rigorous with von Neumann boundary conditions). With normally distributed velocities, the distribution of particles will converge towards a uniform distribution. Further, this model can be trained in a simulation-free manner. We trained this HGF model on a simple toy distribution (see figure 3 and figure 4 in the attached PDF). Such a HGF model is distinct from previous models and illustrates that HGFs are not restricted to only Oscillation HGFs, diffusion models, or flow matching. Let A : R d × R → R d be a time-dependent vector field and ψ the diffeomorphic flow defined by the ODE:d dt x(t) =A(x(t), t) x(0) =x i.e. t → ψ t (x)is a solution to the above ODE. With this, we get:d 2 d 2 t x(t) = d dt d dt x(t) = d dt A(x(t), t) = D x A(x(t), t) ẋ(t) + ∂ ∂t A(x(t), t) = D x A(x(t), t)A(x(t), t) + ∂ ∂t A(x(t), t)Therefore, if we define the force field,F (x, t) = D x A(x, t)A(x, t) + ∂ ∂t A(x, t)we can extend the state space to (x, v) and consider the ODE:(x(0), v(0)) =(x, A(x, 0)) ( d dt x(t), d dt v(t)) =(v(t), F (x(t), t))Then every solution (x t , v t ) to the above ODE is also a solution to the flow matching ODE and vice versa. The conditional velocity predictor loss looks as follows:E[∥V ϕ (ψ t (x t ), t) - d dt ψ t (x t )∥ 2 ]This is exactly the conditional flow matching loss (see equation ( 14) in [31] ).In this section, we discuss the relation between HGF and EDM [26] . The EDM paper assumes the perturbation kernel is isotropic Gaussian with standard deviation σ(t). Thus, the intermediate distribution p t (•; σ(t)) = p data * N (0, σ(t)I). If we further scale the original variable x with s(t) and consider ỹ = s(t)x, [26] shows that the corresponding backward ODE of ỹ is as follows:EQUATIONWe will show that, the minimizer of the objective of the Oscillation HFG, i.e., E y∼π,v∼N (0,I d ) [∥V ϕ (cos(t)y + sin(t)v, t) -[-sin(t)y + cos(t)v]∥ 2 ], equals to the drift term in Eq. 55, when setting s(t) = cos(t) and σ(t) = tan(t).• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 


Decision: This work develops methods for generative modeling based on Hamiltonian dynamics. In this setup the authors combine a learned force field (which approximate the score function of the data distribution) with an auxiliary velocity variable, for which they define a learned Hamiltonian velocity predictor (HVP). The HVP can be used to define a Hamiltonian score matching method, which then leads to the definition of Hamiltonian Generative Flows (HGFs), a class of generative models that recover diffusion models and flow matching as special cases (with zero force field).

Reviewers are in agreement that this is a clearly paper presents an interesting new Hamiltonian perspective on generative modeling. Reviewer do note some weaknesses in experiments, notably the slightly poor performance relative to EDM, though reviewers do not appear to find this a key consideration. Overall, this is a submission that is above the bar for acceptance. In the camera ready, do please make sure to discuss the connections with the AGM method as pointed out by reviewer 6n8r.

Reviews: review: paper_summary: paper_summary: The paper explores the application of Hamiltonian formalism for generative modelling.
In this framework, the score function is interpreted as a force field F, thus optimizing the parameters of F yields a score-matching objective.
Authors build the relation to Flow Matching and Diffusion models as special cases of their framework and suggest another special case, which they validate on an image generation task.
strength_weakness: strength_weakness: - in the general case, the framework requires joint learning of two functions: force field and velocity.
The authors themselves compare the optimization to GANs, which raises the question of training stability.
- the authors do not compare the runtime and memory consumption of the framework with other SOTA methods, e.g. Flow Matching.
- the experimental section is somewhat limited to a toy experiment and a single image-generation task.
It would be highly beneficial to evaluate the framework on another task (e.g. physics-related), although it is not critical given the high theoretical contribution of the work.
- it is not clear how image generation is influenced by Hamiltonian dynamics. What is x, and what is v?
questions: questions: - I am curious if it would be simple to enforce equivariance as in Equivariant Hamiltonian Flows by Rezende et al.
- as the force field form is a choice now, did authors consider alternative forms other than osciallators? I am curious what is the intuition behind the choice.
- authors mentioned conservation properties in the context of physical simulations.
Do you think that they will still be beneficial in the general case, e.g. image/video generation?
Even in physical systems, if the energy is not conserved, would those properties be a limiting factor somehow?
limitations: limitations: Perhaps the main limitation for me is the discussion of the framework's theoretical properties, e.g., conservation of energy and volume, in the context of generative modelling. It would be beneficial to see an analysis of how those inductive biases will play out in specific cases (e.g., physical simulations or image generation). Besides, it would be beneficial to provide computational details (runtime / memory) and compare them to SOTA methods.
paper_summary: paper_summary: In this work, the authors proposed a new generative modeling approach called Hamiltonian Score Matching, which is motivated by classical Hamiltonian mechanics. This approach estimates score functions by augmenting data via Hamiltonian trajectories, and further motivates Hamiltonian Generative Flows. The authors also discuss the design space of force fields in Hamiltonian Generative Flows and connect it with harmonic oscillators. Experiments are conducted to verify the effectiveness of the proposed approach.
strength_weakness: strength_weakness: 1. Regarding the methodology: The authors spend efforts to establish Hamiltonian Generative Flows from Parametrized Hamiltonian ODEs, Hamiltonian velocity predictor, and Theorem 1 for Hamiltonian Score Matching.   Although the benefits such as bringing additional freedoms for force fields and prediction objectives, the scalability of this approach is restricted as stated by the authors. The necessary trajectory simulation and min-max optimization largely limit the usage of the proposed approach. In Section 7, the authors instead seek to use a pre-defined force field instead of parametrized models, which makes the HGFs degenerate to pre-defined ODEs, further indicating the theory-practice gap of the proposed approach.

2. Regarding the experiments: scales and types of experiments in this work are rather limited. In addition to the simulation experiment, the authors only conducted experiments on CIFAR-10, which even shows that the proposed HGFs cannot bring explicit benefits compared to previous approaches. It would largely improve the quality of this work if the authors could provide more diverse and strong experimental results for the proposed approach.

Overall, I think the paper provides an interesting perspective, but the current issues discount the ratings of my evaluation. Currently, I vote for borderline acceptance and I will carefully read the rebuttal and other reviews to decide whether to increase or decrease my scores.
questions: questions: See the above section.
limitations: limitations: The authors discussed the limitations of this work in the Conclusion section.
paper_summary: paper_summary: The authors introduce Hamiltonian velocity predictors (HVPs) as a tool for score matching and generative models. They present two innovations constructed with HVPs: Hamiltonian Score Matching (HSM), a novel generative model that encompasses diffusion models, and flow matching as HGFs with zero force fields. They showcase the extended design space of force fields by introducing Oscillation HGFs, a generative model inspired by harmonic oscillators. They show experiments validating the theoretical insights about HSM.
strength_weakness: strength_weakness: a. Why does it produce inferior performance on image generation compared to EDM?  The author could present the reasons for this. 

b. The authors also need to present the computation cost/efficiency comparison of different methods to validate the proposed method.
questions: questions: See weakness.
limitations: limitations: The authors did not provide a clear discussion of the limitations of the proposed method. This should be included in the paper.
paper_summary: paper_summary: The paper proposes the Hamiltonian Score Matching framework, which is a new general framework of generative models. Inspired by the Hamiltonian dynamics in classical and statistical mechanics, the framework uses the Hamiltonian dynamics to generate data, which is also called Hamiltonian generative flow in the paper. The paper theoretically proves that the force field coincides with the score function if and only if its optimal velocity predictor is zero and proposes the HSM loss function. The paper also shows that the Flow Matching framework and the diffusion model are special cases of the HSM framework. Additionally, a new generative flow called Oscillation Hamiltonian generative flow is constructed. Experiments on several datasets show the experimental performance of the HSM framework.
strength_weakness: strength_weakness: 1. The experimental result is not so good. As the paper writes, the Oscillations HGFs can surpass most baselines but still lack behind the EDM model. I think more experimental results should be provided to show that the HSM framework is a good framework for generative models.
questions: questions: What is the motivation for designing Oscillation HGFs? I'm curious about how to design a force field in the HSM framework.
limitations: limitations: N/A
paper_summary: paper_summary: This paper proposes a new matching method based on Hamiltonian mechanics. The proposed method, algorithmically, is essentially a second-order ODE in which the acceleration (i.e., the drift of the velocity channel) is approximated by deep neural networks. The method is well grounded on the theory of classical Hamiltonian mechanics and statistical physics. As a higher-order system, the proposed method naturally absorbs prior diffusion and flow methods as special cases. Experiments were conducted on a synthetic dataset (Gaussian mixtures) and image generation.
strength_weakness: strength_weakness: - It's unfortunate that, from an algorithmic standpoint, the resulting algorithm---after all these nice theoretical results---remains somewhat similar to a few recently-proposed second-order methods, e.g., AGM (https://arxiv.org/pdf/2310.07805 ICLR'24 oral) or CLD (ref. [12]). 

- The use of harmonic oscillator in Sec 7 is not very motivated in my opinion. There could be other alternatives that also fulfill C1 & C2 conditions. Practically speaking, sampling (x,v) in (23) from sin/cos doesn't seem to encourage straighter generative process.  

- (minor) Typo in L236, I think you meant C2 not C1.

Overall, I have a mixed feeling on this paper --- on one side, I do enjoy, and acknowledge, their theoretical contributions, yet on another side, I feel like the resulting algorithm is rather weak.
questions: questions: - Can the authors comment on the differences compared to AGM? I think AGM also constructs simulation-free (xt, vt) samples then matches the drift of the velocity.

- Have the author ever tried to jointly train both $\theta$ and $\phi$, as mentioned in L179? I understand this may most likely be impractical but just curious from an exploratory standpoint.

- Could the authors elaborate more on L175 the "optimal" choice of T, and the first hypothesis in L184?
limitations: limitations: Limitations were mentioned in Sec 10.

Comments: title: Official Comment by Reviewer g7vb
Author-Reviewer Discussion Phase (ends Aug 13 midnight AOE)
comment: Thank the authors for the response. I will keep the score unchanged.
I thank the authors for the reply. I've decided to keep my score.
Thanks for your reply. All of my concerns have been addressed. I'll keep my positive rating of 6.
I thank the authors for their response and additional experiments, which certainly prove the positive contribution of the submission.
I will retain my score.
Thank you for your clarifications. Most of my concerns have been addressed. I choose to increase my rating to 6.
Dear reviewers,

Thank you for your efforts so far! The authors have responded to each of your reviews and have also provided some additional results in the PDF attached to the overall response. Could those you take a moment to read the author response and follow up regarding any points that may further inform your opinion on the paper?

With thanks, AC