Symplectically Integrated Symbolic Regression
of Hamiltonian Dynamical Systems

Hamiltonian mechanics is a reformulation of classical mechanics introduced in 1833 by William Rowland Hamilton. The Hamiltonian formalism provides a consistent approach for mathematically describing a variety of physical dynamical systems; it is heavily used in quantum mechanics and statistical mechanics, often for systems with many degrees of freedom (Morin, 2008; Reichl, 2016; Taylor, 2005; Sakurai & Napolitano, 2010; Hand & Finch, 1998).

In Hamiltonian mechanics, dynamical systems are governed via their Hamiltonian $\mathcal{H}(\mathbf{p},\mathbf{q})$, a function of the momenta and position vectors ($\mathbf{q}$ and $\mathbf{p}$ respectively) of each object in the system. The value of the Hamiltonian generally measures the total energy of the system; for systems where energy is conserved, $\mathcal{H}$ is time-invariant. Further, we say that $\mathcal{H}$ is separable if it can be decomposed into two additively separable functions $T(\mathbf{p})$ and $V(\mathbf{q})$, where $T$ measures kinetic energy and $V$ measures potential energy.

The time evolution of a physical dynamical system is uniquely defined by the partial derivatives of its Hamiltonian:

The time evolution of a Hamiltonian system is symplectomorphic, preserving the symplectic two-form $d\mathbf{p}\land d\mathbf{q}$. As a result, researchers generally time-evolve Hamiltonians with symplectic integrators, which are designed to preserve this two-form. We employ a fourth-order symplectic integrator (pseudocode in Appendix A) (Forest & Ruth, 1990).

Hamiltonians commonly exhibit coupling: they can be written as a series of additively separable functions where each function takes some smaller subset of variables. Consider some arbitrary n-body, m-dimensional separable Hamiltonian $\mathcal{H}(\mathbf{p},\mathbf{q})=T(\mathbf{p})+V(\mathbf{q})$ where the bodies are labeled $1,\dots,n$ and the dimensions are labeled $1,\dots,m$. What coupling might occur?

Without loss of generality, consider $T(\mathbf{p})$. Denote the momentum of body $i$ in dimension $j$ as $p_{i,j}$. We say that $T(\mathbf{p})$ exhibits complete decoupling if it may be rewritten as the sum of functions that each take only one variable in the system, i.e. it may be written in the form:

We say that $T(\mathbf{p})$ exhibits dimensional coupling if it may be rewritten as the sum of functions that each take all dimensional variables for only a single particle in the system, i.e. it may be written in the form:

Finally, we say that $T(\mathbf{p})$ exhibits pairwise coupling if it may be rewritten as the sum of functions that each take all dimensional variables for a pair of particles in the system, i.e. where $\mathbf{p_{i}}$ denotes all momenta variables of particle $i$, $T$ may be written in the form:

Consider some arbitrary coupling function $f(x_{1},\dots,x_{n})$. Often, $f$ may be written as a composition of functions $f(x_{1},\dots,x_{n})=g(h(x_{1},\dots,x_{n}))$ where $h$ performs one of several common reduction operations, such as summing, multiplication, Euclidean distance, or Manhattan distance. In these instances, we would say that particles are coupled via the given inner operation, such as their sum. If $h$ is known, then one may learn $f$ by learning $g$, which tends to be much simpler than learning the entire $f$ directly.

Suppose that a Hamiltonian function exhibits pairwise coupling, i.e., it may be written as the sum of many additively separable coupling functions. Often, across all pairs of particles, these coupling functions are identical (symmetric). For instance, in an n-body gravitational system with particles of equal mass, or an n-body electrostatic system with particles of equal charge, all particles are pairwise coupled using the exact same coupling function.

As detailed in Sections 3.2, 3.3, and 3.4, Hamiltonians commonly exhibit coupling properties. These properties may be used to simplify the Hamiltonian discovery process by imposing restrictions that dramatically reduce the solution space. Thus, the first step of SISR is to find what coupling is taking place in the Hamiltonian, whether it can be reduced into simpler composite functions, and whether it is symmetric across particles or pairs of particles.

Our approach for doing this is rooted in symplectic neural networks as effective black-box predictors of Hamiltonian dynamical systems. Specifically, we can construct multi-layer neural networks that enforce specific coupling properties. Then, we may train these networks by using them to time-evolve the dynamical system through a symplectic integrator. By seeing how these networks perform with and without imposed coupling properties, we may deduce whether those properties hold for the system. To the best of our knowledge, we are the first work to extract coupling properties in the context of Hamiltonian dynamical systems, and we are not aware of any other work that employs symplectic neural networks in a similar manner.

With our coupling properties determined, we are ready to begin generating candidate symbolic Hamiltonians. We do this via a deep symbolic regression approach, employing an LSTM-RNN.

First, note that free-form symbolic functions may be represented as trees where internal nodes are functional operators and leaves are variables or constants. The pre-order traversal of this expression tree produces a sequence that uniquely defines the symbolic expression. Thus, sequence-generating RNNs can be used in an autoregressive manner to produce symbolic functions. We refer to functional operators, variables, and constants as tokens in the sequence.

To generate a symbolic expression, an LSTM-RNN generates a sequence from left-to-right. As input, the network takes a one-hot encoding of the parent and sibling token for the current token being sampled. It then emits a softmax distribution, which is adjusted to eliminate all invalid tokens. This distribution may be sampled to generate the next token in the sequence. We also employ an epsilon-greedy mutation approach where the next token is sampled uniformly from the list of valid tokens with $\mu$ probability, where $\mu$ is a pre-defined hyperparameter. This process is depicted in Step 1 of Figure 2.

Due to our use of a fourth-order symplectic integrator, we require that any Hamiltonian generated by the LSTM-RNN is of the form $\mathcal{H}(\mathbf{p},\mathbf{q})=T(\mathbf{p})+V(\mathbf{q})$. This constraint is applied in situ, meaning that the softmax distribution over the operators is constantly adjusted as expressions are generated so that only Hamiltonians of this form will ever be produced. This same approach is used to incorporate any coupling properties previously discovered: the emitted distribution is adjusted prior to sampling a token so that we will only ever see Hamiltonians that exhibit the desired coupling. To the best of our knowledge, this is the first instance of an LSTM-RNN being used to generate separable Hamiltonians with desired coupling properties.

Once a Hamiltonian sequence is generated, it is translated into an executable model with auto-differentiation (Steps 2 and 3 of Figure 2). The first reason for this is to that constants can be optimized–the sequence contains placeholder “ephemeral” constants that do not refer to specific values. The second purpose is so that the performance of the expression can be assessed. This is done by numerically integrating the Hamiltonian to predict the system, which requires the ability to take its gradients.

As mentioned previously, candidate Hamiltonian functions must have their placeholder constants trained. This is performed via a symplectic optimization approach (Step 4 of Figure 2). Using a fourth-order symplectic integrator and a time-series of observation data for the system, we use the candidate Hamiltonian to time-evolve the system and obtain predicted future states. We then take the loss between our predicted future states and the actual future states contained in the observation data, backpropagating through our integrator to train the constants and thus imparting an energy-conserving symplectic prior onto any Hamiltonian generated by the RNN. This approach, which has been applied to black-box system prediction methods with great success, has yet to be employed for any free-form symbolic regression model (Chen et al., 2019; DiPietro et al., 2020).

Thus far, we have presented an RNN that can generate separable Hamiltonians (with specific coupling properties) that are then trained via backpropagation through a symplectic integrator. However, in order for a symbolic Hamiltonian search to succeed, the RNN must be trained to produce progressively better-fitting Hamiltonians.

To do this, we employ the risk-seeking policy gradients approach presented by Petersen et al. (2019), which maximizes the best-case expectation of rewards for expressions generated by the policy (in this case the RNN). To do this, we train on the top $\epsilon$ quantile in some given batch, where $1-\epsilon$ is referred to as the risk factor. First, let this quantile consist of $M$ expressions where $\tau_{i}$ denotes the $i$-th expression in the quantile. If $f$ defines a symplectic integration function, $R$ defines a reward function, $R_{\epsilon}$ defines the reward of the top $\epsilon$ quantile, and $J(\theta)$ defines the expectation of rewards in the top $\epsilon$ quantile conditioned on the weights $\theta$ of the RNN, we have:

In the event that the expression generation involved mutation, we compute $\log p(\tau_{i}|\theta)$ using the probabilities produced by the RNN–not the uniform probabilities generated to sample the mutated operator.

Further, as in Petersen et al. (2019), we adopt the maximum entropy framework and include an entropy term in the loss (Haarnoja et al., 2018). Where $\lambda_{H}$ is the entropy coefficient, we have:

Note that $H$ refers to the Shannon entropy function and is not a Hamiltonian.

This training process corresponds to Step 5 of Figure 2. We are not aware of any other approach for Hamiltonian dynamical systems that employs these reinforcement learning-oriented training approaches.

The harmonic oscillator is a one-dimensional, single-particle dynamical system described by the Hamiltonian

where $m$ is the mass of the particle and $\omega$ refers to its oscillator frequency.

We trained SISR on 3 different datasets of harmonic oscillator motion with varying initial conditions, each consisting of only $30$ points from $t=0$ to $t=3$. In addition to variables and constants, SISR had access to the operators $\{+,-,\div,\cdot,\wedge\}$. The $T$ and $V$ of each generated Hamiltonian were limited to be between $1$ and $8$ operators long. SISR was applied to each dataset 5 times, recovering the correct Hamiltonian for all 15 trials. This took an average of $4.60\pm 3.25$ batches, or $225.83\pm 101.99$ seconds. In 4 out of 15 instances, the Hamiltonian was uncovered after a single epoch. Next, we applied SISR to the same datasets with $(\mu=0,\sigma=0.001)$ of Gaussian noise. Once again, it uncovered the correct Hamiltonian in all trials, requiring an average of $5.00\pm 3.34$ batches, or $252.76\pm 109.04$ seconds.

The pendulum is an angular, single-particle dynamical system described by the Hamiltonian

where $m$ is mass, $g$ is gravity, and $\ell$ is length.

Using the functional operators $\{+,-,\div,\cdot,\wedge,\cos,\sin\}$, we trained SISR on 3 different datasets of pendulum motion with varying initial conditions, containing $50$ points from $t=0$ to $t=5$. The $T$ and $V$ of each generated Hamiltonian were limited to be between $1$ and $8$ operators long. With no noise present, SISR recovered the correct Hamiltonian in all 15 trials, taking an average of $21.80\pm 11.18$ batches, or $1049.41\pm 465.10$ seconds. When $(\mu=0,\sigma=0.005)$ of Gaussian noise was added to each dataset, SISR recovered all correct Hamiltonians in $28.73\pm 13.68$ batches, or $1117.95\pm 497.57$ seconds. Hamiltonians generated by SISR for the pendulum problem are depicted in Figure 3.

The two-body gravitational system is a two-dimensional dynamical system described by the Hamiltonian

with mass $m$ for particles $i$ and $j$. Note that we work explicitly with two-body systems where all masses are equal and $G=1$.

We began by extracting coupling properties via the symplectic neural network procedures described in Section 4.1. Using this process, we first determined that particle positions were pairwise coupled, as well as that particle momenta were completely decoupled. Then, we were able to decompose the position coupling into a composite function of Euclidean distance. Finally, we determined that the functions acting on the momenta of each particle are symmetric. The exact numerical results used to guide this search process are available in Appendix B.3. Once these coupling constraints are imposed on the RNN, the solution space shrinks dramatically. Rather than learning $V(\mathbf{q})$ and $T(\mathbf{p})$ outright as in our previous experiments, we are effectively learning the pairwise position coupling function between particles, which is itself a function of their Euclidean distance, and the decoupled momentum function acting on each particle.

Using the functional operators $\{+,-,\div,\cdot,\wedge\}$, we trained SISR on 3 different datasets of two-body motion with varying initial conditions, containing $200$ points from $t=0$ to $t=20$. The $T$ and $V$ of each generated Hamiltonian were limited to be between $1$ and $12$ operators long. With no noise present, SISR recovered the correct Hamiltonian for all 15 trials, taking an average of $1.47\pm 1.30$ batches, or $227.42\pm 62.94$ seconds. When $(\mu=0,\sigma=0.001)$ of Gaussian noise was added to each dataset, SISR recovered all correct Hamiltonians in $1.80\pm 1.15$ batches, or $236.79\pm 55.45$ seconds. Note how powerful our extracted coupling priors are: with them imposed, the two-body Hamiltonian is learned more quickly than a pendulum Hamiltonian.

The three-body gravitational system is a two-dimensional dynamical system described by the Hamiltonian

with mass $m$ for particles $i$, $j$, and $k$. This system exhibits chaotic motion, meaning that small perturbations in initial conditions diverge exponentially in the future (Sprott, 2010).

We followed a nearly identical coupling extraction process to our two-body experiments, identifying that all particles are pairwise coupled as a function of their Euclidean distance, as well as all momenta being completely decoupled and symmetric. These results are detailed in Appendix B.3.

Using the functional operators $\{+,-,\div,\cdot,\wedge\}$, we trained SISR on 3 different datasets of three-body motion with $200$ points each. The $T$ and $V$ of each generated Hamiltonian were limited to be between $1$ and $18$ operators long. Five trials were conducted on each dataset. With no noise present, SISR took an average of $2.40\pm 2.32$ epochs, or $588.92\pm 238.49$ seconds, to recover the Hamiltonian (all were recovered). With $(\mu=0,\sigma=0.001)$ of Gaussian noise present, SISR required an average of $2.13\pm 1.77$ epochs, or $575.48\pm 201.20$ seconds (all were recovered).

We benchmarked SISR against two other methods for obtaining physical governing equations for data: SINDy and SSINNs. Both of these methods rely upon sparse regressions in user-defined function spaces. They obtained reasonable NRMSE for many problems but only extracted correct governing equations in 2 of 8 instances. Indeed, with so little data, it is easy for sparse-regression oriented methods to overfit to excess terms. Finally, we compared to the vanilla deep symbolic regression model, which we used to generate systems of differential equations instead of Hamiltonians. These results are summarized in Table 1 with further experimental details available in Appendix B.4.