Hamiltonian Mechanics ¶

Hamiltonian mechanics is a way to describe how physical systems, like planets or pendulums, move over time, focusing on energy rather than just forces. By reframing complex dynamics through energy lenses, this 19th-century physics framework now powers cutting-edge generative AI. It uses generalized coordinates \( q \) (like position) and their conjugate momenta \( p \) (related to momentum), forming a phase space that captures the system's state. This approach is particularly useful for complex systems with many parts, making it easier to find patterns and conservation laws.

Mathematical Reformation: From Second-Order to Phase Flow ¶

Newton's \( F = m\ddot{q} \) requires solving second-order differential equations, which become unwieldy for constrained systems or when identifying conserved quantities.

Lagrangian Prelude: Action Principles¶

The Lagrangian \( \mathcal{L}(q, \dot{q}) = K - U \) leads to Euler-Lagrange equations via variational calculus: $$ \frac{d}{dt}\left( \frac{\partial \mathcal{L}}{\partial \dot{q}} \right) - \frac{\partial \mathcal{L}}{\partial q} = 0 $$

But these remain second-order. The critical leap comes through Legendre Transformation \( (\dot{q} \rightarrow p) \). The Hamiltonian is derived from the Lagrangian through a Legendre transformation by defining the conjugate momentum as \( p_i = \frac{\partial \mathcal{L}}{\partial \dot{q}_i} \); then the Hamiltonian can be written as: $$ H(q,p) = \sum_i p_i \dot{q}_i - \mathcal{L}(q, \dot{q}) $$

We can write \( H(q,p) \) more intuitively as: $$ H(q,p) = K(p) + U(q) $$

This flips the script: instead of \( \dot{q} \)-centric dynamics, we get symplectic phase flow.

This transformation is not canonical because it does not preserve the Poisson bracket structure.

Newton vs. Lagrange vs. Hamilton: A Philosophical Showdown¶

Hamilton's Equations: The Geometry of Phase Space ¶

The phase space is a mathematical space where we can represent the set of possible states of a physical system. For a system with \( n \) degrees of freedom, the phase space is a \( 2n \)-dimensional space, often visualized as a map where each point \( (q, p) \) represents a unique state. The evolution of the system is described by the motion of a point in this space, governed by Hamilton's equations.

This formulation offers several advantages. It makes it straightforward to identify conserved quantities and symmetries through canonical transformations and Poisson brackets, which provides deeper insights into the system's behavior. For instance, Liouville's theorem states that the volume in phase space occupied by an ensemble of systems remains constant over time, expressed as:

or equivalently:

where \( \rho(q, p, t) \) is the density function. This helps us to represent the phase space flows and how they preserve area under symplectic transformations. Its relation to symplectic geometry enables mathematical properties that are directly relevant to many numerical methods. For instance, it enables Hamiltonian Monte Carlo to perform well in high-dimensions by defining MCMC strategies that increases the chances of accepting a sample (particle).

Symplecticity: The Sacred Invariant¶

Hamiltonian flows preserve the symplectic 2-form \( \omega = \sum_i dq_i \wedge dp_i \).

import torch

# Conceptual representation of dq and dp as tensors
dq = torch.tensor([1.0])  
dp = torch.tensor([1.0])  

# "Wedge product" conceptually represented using an outer product
omega = torch.outer(dq, dp) - torch.outer(dp, dq)

print(omega)

Numerically, this means good integrators must respect:

How are We Going to Solve it?¶

Numerical solvers for differential equations inevitably introduce errors that affect solution accuracy. These errors manifest as deviations from the true trajectory in phase space, particularly noticeable in energy-conserving systems like the harmonic oscillator. The errors fall into two main categories: local truncation error, arising from the approximation of continuous derivatives with discrete steps (proportional to \( \mathcal{O}(\epsilon^n+1) \) where \( \epsilon \) is the step size and n depends on the method); and global accumulation error, which compounds over integration time.

To overcome this, we turn to symplectic integrators—methods that respect the underlying geometry of Hamiltonian systems, leading us naturally to the Leapfrog Verlet method, a powerful symplectic alternative. 🚀

Symplectic Numerical Integrators ¶

Leapfrog Verlet¶

For a separable Hamiltonian \( H(q,p) = K(p) + U(q) \), where the corresponding probability distribution is given by:

the Leapfrog Verlet integrator proceeds as follows:

This Störmer-Verlet scheme preserves symplecticity exactly, with local error \( \mathcal{O}(\epsilon^3) \) and global error \( \mathcal{O}(\epsilon^2) \). You can read more about numerical methods and analysis in Python here.

Why Symplectic Matters¶

They're the reversible neural nets of physics simulations!

Symplectic integrators like Leapfrog Verlet are critical for long-term stability in Hamiltonian systems.

• Phase space preservation: The volume in \( (q, p) \)-space is conserved exactly, avoiding artificial energy drift.
• Approximate energy conservation: While energy \( H(q,p) \) is not perfectly conserved (due to \( \mathcal{O}(\epsilon^2) \) error), it oscillates near the true value over exponentially long timescales.
• Practical relevance: This makes symplectic integrators indispensable in molecular dynamics and Hamiltonian Monte Carlo (HMC), where accurate sampling relies on stable trajectories.

Euler's method (first-order) systematically injects energy into the system, causing the characteristic outward spiral seen in the plots. Modified Euler's method (second-order) significantly reduces this energy drift. Most importantly, symplectic integrators like the Leapfrog method preserve the geometric structure of Hamiltonian systems even with relatively large step sizes by maintaining phase space volume conservation. This structural preservation is why Leapfrog remains the preferred method for long-time simulations in molecular dynamics and astronomy, where energy conservation is critical despite the visible polygon-like discretization artifacts at large step sizes.

Non-symplectic methods (e.g., Euler-Maruyama) often fail catastrophically in these settings.

Integrator	Symplecticity	Order	Type	Local Error	Global Error	Suitable For	Computational Cost
Euler Method		1	Explicit	O(ε²)	O(ε)	Quick prototypes and Short-term simulations of general ODEs	Low
Symplectically Euler		1	Explicit	O(ε²)	O(ε)	Simple Hamiltonian systems	Low
Leapfrog (Verlet)		2	Explicit	O(ε³)	O(ε²)	Molecular dynamics, Long-term simulations of Hamiltonian systems	Moderate
Runge-Kutta 4		4	Explicit	O(ε⁵)	O(ε⁴)	Short-term accuracy, General ODEs, but not recommended for long-term Hamiltonian systems	High
Forest-Ruth Integrator		4	Explicit	O(ε⁵)	O(ε⁴)	High-accuracy long-term simulations	High
Yoshida 6th-order		6	Explicit	O(ε⁷)	O(ε⁶)	High-accuracy	High
Heun’s Method (RK2)		2	Explicit	O(ε³)	O(ε²)	General ODEs requiring moderate accuracy	Moderate
Third-order Runge-Kutta		3	Explicit	O(ε⁴)	O(ε³)	When higher accuracy than RK2 is needed without the cost of RK4	High
Implicit Midpoint Rule		2	Implicit (solving equations)	O(ε³)	O(ε²)	Hamiltonian systems, stiff problems	High
Fourth-order Adams-Bashforth		4	Multi-step (explicit)	O(ε⁵)	O(ε⁴)	Non-stiff problems with smooth solutions, after initial steps	Low
Backward Euler Method		1	Implicit (solving equations)	O(ε²)	O(ε)	Stiff problems, where stability is crucial	High

Hamiltonian Monte Carlo¶

Hamiltonian Monte Carlo (HMC) is a Markov chain Monte Carlo (MCMC) method that leverages Hamiltonian dynamics to efficiently sample from complex probability distributions, particularly in Bayesian statistics and machine learning.

From Phase Space to Probability Space¶

HMC interprets target distribution \( P(z) \) as a Boltzmann distribution:

Substituting into this formulation, the Hamiltonian gives us a joint density:

where \( p(q|D) \) is the likelihood of the given data \( D \) and T=1 and therefore removed. We estimate our posterior distribution using the potential energy \( U(q) \) since \( P(q,p) \) consists of two independent probability distributions.

Augment with artificial momentum \( p \sim \mathcal{N}(0,M) \), then simulate Hamiltonian dynamics to propose new \( q' \) based on the distribution of the position variables \( U(q) \) which acts as the "potential energy" of the target distribution \( P(q) \), thereby creating valleys at high-probability regions.

For more on HMC, check out this explanation or this tutorial.

The kinetic energy with the popular form of \( K(p) = \frac{1}{2}p^T M^{-1} p \), often Gaussian, injects momentum to traverse these landscapes. Crucially, the mass matrix \( M \) plays the role of a preconditioner - diagonal \( M \) adapts to parameter scales, while dense \( M \) can align with correlation structure. \( M \) is symmetric, positive definite and typically diagonal.

What is Positive Definite?

Positive Definite: For any non-zero vector \( x \), the expression \( x^T M x \) is always positive. This ensures stability and efficiency.

\[ x^T M x > 0 \]

The HMC Algorithm¶

The algorithm involves:

1. Initialization: Start with an initial position \( q_0 \) and sample momentum \( p_0 \sim \mathcal{N}(0,M) \).

2. Leapfrog Integration: Use the leapfrog method to approximate Hamiltonian dynamics. For a step size \( \epsilon \) and L steps, update:

3. Half-step momentum: \( p(t + \frac{\epsilon}{2}) = p(t) - \frac{\epsilon}{2} \frac{\partial U}{\partial q}(q(t)) \)

4. Full-step position: \( q(t + \epsilon) = q(t) + \epsilon \frac{\partial K}{\partial p}(p(t + \frac{\epsilon}{2})) \), where \( K(p) = \frac{1}{2} p^T M^{-1} p \), so \( \frac{\partial K}{\partial p} = M^{-1} p \)
5. Full-step momentum: \( p(t + \epsilon) = p(t + \frac{\epsilon}{2}) - \frac{\epsilon}{2} \frac{\partial U}{\partial q}(q(t + \epsilon)) \)

This is repeated L times to get proposed \( \dot{q} \) and \( \dot{p} \).

1. Metropolis-Hastings Acceptance: Accept the proposed \( \dot{q} \) with probability \( \min(1, e^{H(q_0,p_0) - H(\dot{q},\dot{p})}) \), where \( H(q,p) = U(q) + K(p) \).

This process generates a Markov chain with stationary distribution \( P(q) \), leveraging Hamiltonian dynamics to take larger, more efficient steps compared to random-walk methods.

This coupled system drives \( (q,p) \) along iso-probability contours of \( P(q) \), with momentum rotating rather than resetting at each step like in Random Walk Metropolis--think of following mountain paths instead of wandering randomly! The key parameters - integration time \( \tau = L\epsilon \) and step size \( \epsilon \) - balance exploration vs. computational cost:

• Short \( \tau \): Local exploration, higher acceptance
• Long \( \tau \): Global moves, risk of U-turns (periodic orbits)

Connection with Energy-Based Models¶

Energy-based models (EBMs) are a class of generative models that define a probability distribution over data points using an energy function. The probability of a data point is proportional to \( e^{-E(x)} \), where \( E(x) \) is the energy function. This formulation is directly analogous to the Boltzmann distribution in statistical physics, where the probability is related to the energy of a state. In Hamiltonian mechanics, the Hamiltonian function \( H(q, p) \) represents the total energy of the system, and the probability distribution in phase space is given by \( e^{-H(q,p)/T} \), where \( T \) is the temperature.

In EBMs, Hamiltonian Monte Carlo (HMC) is often used to sample from the model's distribution. HMC leverages Hamiltonian dynamics to propose new states, which are then accepted or rejected based on the Metropolis-Hastings criterion. This method is particularly effective for high-dimensional problems, as it reduces the correlation between samples and allows for more efficient exploration of the state space. For instance, in image generation tasks, HMC can sample from the distribution defined by the energy function, facilitating the generation of high-quality images.

EBMs define probability through Hamiltonians:

Potential Research Directions¶

Symplecticity in Machine Learning Models¶

Incorporate the symplectic structure of Hamiltonian mechanics into machine learning models to preserve properties like energy conservation, which is crucial for long-term predictions. Generalizing Hamiltonian Neural Networks (HNNs), as discussed in Hamiltonian Neural Networks, to more complex systems or developing new architectures that preserve symplecticity

HMC for Complex Distributions:¶

HMC for sampling from complex, high-dimensional, and multimodal distributions, such as those encountered in deep learning. Combining HMC with other techniques, like parallel tempering, could handle distributions with multiple modes more effectively.

Combining Hamiltonian Mechanics with Other ML Techniques:¶

Integrate Hamiltonian mechanics with reinforcement learning to guide exploration in continuous state and action spaces. Using it to model the environment could improve exploration strategies, as seen in potential applications in robotics. Additionally, using Hamiltonian mechanics to define approximate posteriors in variational inference could lead to more flexible and accurate approximations.

Hamiltonian GANs¶

Employing Hamiltonian formalism as an inductive bias for the generation of physically plausible videos with neural networks. Imagine generator-discriminator dynamics governed by:

The possibilities make phase space feel infinite...

Useful Links¶

• An Introduction to Multistep Methods: Leap-frog
• The beginners guide to Hamiltonian Monte Carlo
• Hamiltonian Monte Carlo
• Hamiltonian Monte Carlo - Stan - Stan explained
• Hamiltonian Mechanics For Dummies: An Intuitive Introduction.
• Hamiltonian mechanics Wikipedia page
• An introduction to Lagrangian and Hamiltonian mechanics Lecture notes

• Hamiltonian Mechanics - Jeremy Tatum, University of Victoria

• Hamiltonian Neural Networks - Blog

• Hamiltonian Neural Networks
• Other: Natural Intelligence - A blog by Sam Greydanus - Many interesting topics