Learning and Interpreting Potentials for Classical Hamiltonian Systems
Harish S. Bhat
cs.LG math.DS physics.class-ph stat.ML










H. S. Bhat



Applied Mathematics Unit, University of California, Merced, USA
hbhat@ucmerced.edu


Learning and Interpreting Potentials for Classical Hamiltonian Systems
    Harish S. Bhat 0000-0001-7631-1831 This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and was supported by the LLNL-LDRD Program under Project No. 19-ERD-009. LLNL-PROC-779792.
    December 30, 2023
========================================================================================================================================================================================================================================================================================




We consider the problem of learning an interpretable potential energy function from a Hamiltonian system's trajectories.  We address this problem for classical, separable Hamiltonian systems.  Our approach first constructs a neural network model of the potential and then applies an equation discovery technique to extract from the neural potential a closed-form algebraic expression.  We demonstrate this approach for several systems, including oscillators, a central force problem, and a problem of two charged particles in a classical Coulomb potential.    Through these test problems, we show close agreement between learned neural potentials, the interpreted potentials we obtain after training, and the ground truth.  In particular, for the central force problem, we show that our approach learns the correct effective potential, a reduced-order model of the system.








§ INTRODUCTION

As a cornerstone of classical physics, Hamiltonian systems arise in numerous settings in engineering and the physical sciences.  Common examples include coupled oscillators, systems of particles/masses subject to classical electrostatic or gravitational forces, and rigid bodies.   For integer d ≥ 1, let (t) ∈ℝ^d and (t) ∈ℝ^d denote, respectively, the position and momentum of the system at time t.  Let T and V denote kinetic and potential energy, respectively.  Our focus here is on classical, separable systems that arise from the Hamiltonian

    H(, ) =T() + V().

In this paper, we consider the problem of learning or identifying the potential energy V from data ( (t), (t) ) measured at a discrete set of times.  We assume T is known.  To motivate this problem, consider the setting of m interacting particles in three-dimensional space; here d = 3m.  Suppose that we are truly interested in a reduced set of variables, e.g., the position and momentum of one of the m particles.  Let us denote the reduced-order quantities of interest by ((t), (t)) ∈ℝ^2 d.   The direct approach is to integrate numerically the 6m-dimensional system of differential equations for the full Hamiltonian (<ref>) and then use the full solution ( (t), (t) ) to compute  ((t), (t)).  While such an approach yields numerical answers, typically, it does not explain how the reduced-order system evolves dynamically in time.  If we suspect that ((t), (t)) itself satisfies a Hamiltonian system, we can search for a potential V() that yields an accurate, reduced-order model for ((t), (t)). If V is interpretable, we can use it to explain the reduced system's dynamics—here we mean interpretability in the sense of traditional models in the physical sciences, which are written as algebraic expressions, not as numerical algorithms.  We can also use the reduced-order model to simulate ((t), (t)) directly, with computational savings that depend on d/d.

There is a rapidly growing literature on machine learning of potential energies in computational/physical chemistry, e.g., <cit.>.  As in these studies, the present work uses neural networks to parameterize the unknown potential.  A key difference is that, in the present work, we apply additional methods to interpret the learned neural potential.  There exists a burgeoning, recent literature on learning interpretable dynamical systems from time series, e.g., <cit.>, to cite but a few.  We repurpose one such method—SINDy (sparse identification of nonlinear dynamics)—to convert the learned neural potential into a closed-form algebraic expression that is as interpretable as classical models.  We apply only one such method for accomplishing this conversion into an algebraic expression; we hope that the results described here lead to further investigation in this area.



§ APPROACH

Assume T() = ∑_i=1^d M_ii^-1 p_i^2 where M is a diagonal mass matrix.  Then, from (<ref>), we can write Hamilton's equations:



    = M^-1
       = -∇ V().


Let the training data consist of a set of R trajectories; the j-th such trajectory is {^j_i, ^j_i }_i=0^N.  Here (^j_i, ^j_i)   denotes a measurement of ((t), (t)) at time t = i h for fixed h > 0.  We choose this equispaced temporal grid for simplicity; this choice is not essential.    Because we treat the kinetic energy T as known, we assume that the training data consists of (possibly noisy) measurements of a system that satisfies (<ref>).  We now posit a model for V that depends on a set of parameters .  For instance, if we model V using a neural network,  stands for the collection of all network weights and biases.  Then we use (<ref>) to  form an empirical risk loss 

    L() = 1/R N∑_j=1^R∑_i=0^N-1^j_i+1 - ^j_i/h + ∇_ V(^j_i; ) ^2.

Let τ = Nh denote the final time in our grid.  Note that (<ref>) approximates E [ (1/τ) ∫_t=0^t=τ(t) + ∇_ V((t); ) ^2 dt], the expected mean-squared error of a random trajectory ((t),(t)) assumed to satisfy (<ref>).

We model V using a dense, feedforward neural network with L ≥ 2 layers.  Because we train with multiple trajectories, the input layer takes data in the form of two tensors—one for  and one for —with dimensions N × R × d.  The network then transposes and flattens the data to be of dimension N R × d.  Thus begins the potential energy function part of the network (referred to in what follows as the neural potential), which takes a d-dimensional vector as input and produces a scalar as output.  Between the neural potential's d-unit input layer and 1-unit output layer, we have a number of hidden layers.  In this model, we typically choose hidden layers to all have ν units where 1 ≤ν≤ d.  As these architectural details differ by example, we give them below.

Note that the loss (<ref>) involves the gradient of V with respect to the input 𝐪.  We use automatic differentiation to compute this gradient.  More specifically, in our IPython/Jupyter notebooks (linked below), we use the  method in TensorFlow.  This  is easy to implement, fast, and accurate up to machine precision.

The trained network gives us a neural potential V : ℝ^d →ℝ.  To interpret V, we apply the SINDy method <cit.>.  We now offer a capsule summary of this technique.  Suppose we have a grid {𝐱^k}_k=1^K of points in ℝ^d.  We use the notation 𝐱^k = (x_1^k, …, x_d^k).  We evaluate V on the grid, resulting in a vector of values that we denote by 𝐕.  We also evaluate on the grid a library of J candidate functions ξ^j : ℝ^d →ℝ for 1 ≤ j ≤ J; each such evaluation results in a vector Ξ^j that we take to be the j-th column of a matrix Ξ.  In d=1, examples of candidate functions are {1, x, x^2, …} or {1, x^-1, x^-2, …}.  In d=2, an example is {1, x_1, x_2, x_1^2, x_1 x_2, x_2^2, …}.  Each candidate function is simply a scalar-valued function on ℝ^d.

Equipped with the K × 1 vector 𝐕 and the K × J matrix Ξ, we solve the regression problem

    𝐕 = Ξ + ϵ

for the J × 1 vector  using an iteratively thresholded least-squares algorithm.  The algorithm has one constant hyperparameter, λ > 0.  The algorithm is then succinctly described as follows: (i) estimate  using ordinary least squares, and then (ii) reset to zero all components of  that are less than the threshold λ.  Once components of  are reset to zero, they stay frozen at zero.  We then repeat steps (i) and (ii) until  stabilizes to its converged value.

As shown recently  <cit.>, this algorithm converges in a finite number of steps to an approximate minimizer of 𝐕 - Ξ^2 + λ^2 _0.  Here _0 denotes the number of nonzero entries of .  Hence, increasing the parameter λ leads to a more sparse set of coefficients .  Once we have fit the regression model in this way, we obtain an interpretable model of V, specifically:

    V(𝐱) = ∑_j=1^Jβ_j ξ^j(𝐱) + ϵ.

If  is highly sparse, most of the coefficients β_j will be zero.  Suppose that the candidate functions ξ^j are well-known functions such as positive or negative powers of the coordinates x_i of the input 𝐱.  In this case, the right-hand side will be a relatively short algebraic expression that is as interpretable as most potential energy functions routinely encountered in classical physics.  The norm of ϵ here captures the error in this sparse approximation of V.  In general, one chooses λ to balance the sparsity of  with the quality of the approximation ϵ.


§ TESTS

We now describe a series of increasingly complex tests that demonstrate the proposed method.  For each such model, we use either exact solutions or fine-scale numerical integration to create a corpus of time series measurements.  Using the time series, we train a neural potential energy model, which we then interpret using SINDy.  We use NumPy/SciPy or Mathematica for all data generation, TensorFlow for all neural network model development/training, and the  package <cit.> in R to interpret the neural potential.  In what follows, the mass matrix M in (<ref>) is the identity unless specified otherwise.  In all cases, we train the neural network using gradient descent.  We are committed to releasing all code/data at <https://github.com/hbhat4000/learningpotentials>.



 §.§ Simple Harmonic Oscillator

The first model we consider is the simple harmonic oscillator (d=1) with Hamiltonian

    H(q,p) = p^2/2 + q^2/2.

Exact trajectories consists of circles centered at the origin in (q,p) space.  For training data, we use R=10 such circles; for 1 ≤ i ≤ 10, the i-th circle passes through an initial condition (q(0),p(0)) = (0, i).  We include N=1000 steps of each trajectory, recorded with a time step of 0.01, in the training set.  Here our goal is to check how closely the neural potential V(q) can track the true potential V(q) = q^2/2.  We take the neural potential model to have two hidden layers, each with 16 units and tanh activations.  We train for 50000 steps at a learning rate of 0.01.



In Figure <ref>, we plot both the trained neural potential V (in red) and the true potential V(q) = q^2/2 (in black).  When plotting V, we have subtracted a constant bias (the minimum obtained value of V) so that the curve reaches a minimum value of zero.  Note that this constant bias is completely unimportant for physics, as only ∇ V appears in Hamilton's equations (<ref>) and the loss (<ref>).  However, the constant bias does lead us to include an intercept in the regression model (<ref>), i.e., to include 1 in the set of candidate functions for SINDy.  We follow this practice in all uses of SINDy below.

We apply SINDy to the learned potential V(q) with candidate functions {1, q, q^2, q^3}.  In the following test, and in fact throughout this paper, we start with λ = 1 and tune λ downward until the error ϵ—between the neural network potential  V(q) and the SINDy-computed approximation—drops below 10^-10.  We find that with λ = 0.04, the estimated system is

    V(q) ≈β_0 + β_2 q^2

with β_0 ≈ -49.18 and β_2 ≈ 0.4978. We see that V closely tracks the true potential V(q) = q^2/2 up to the constant bias term, which can be ignored.







 §.§ Double Well

Let us consider a particle in a double well potential (d=1)

    V(q) = x^2 (x-2)^2 - (x-1)^2.

We take the kinetic energy to be T(p) = p^2/2.    We now use explicit Runge-Kutta integration in Mathematica to form three training sets:


  * Training set 𝒯_1 includes R=10 trajectories with random initial conditions (q(0), p(0)) chosen uniformly on [-1,1]^2, one of which has sufficiently high energy to visit both wells.

  * Training set 𝒯_2 consists of R=2 trajectories, each of which starts and stays in an opposing well.  The first trajectory has initial condition (q(0),p(0)) = (3,0) while the second has initial condition (q(0),p(0)) = (-1,0).  These q(0) values are symmetric across q=1, the symmetry axis of V(q).

  * Training set 𝒯_3 has only R=2 trajectories that stay in the left well only.

For each trajectory, we record 5001 points at a time step of 0.001.  We take the neural potential model to have two hidden layers, each with 16 units and tanh activations.  For each training set 𝒯_m, we train for 50000 steps at a learning rate of 0.01.

We seek to understand how the choice of training set 𝒯_m affects the ability of the neural potential V to track the true potential (<ref>). We plot and discuss the results in Figures <ref> and <ref>.  Overall, the neural potentials trained using 𝒯_1 and 𝒯_2 match V(q) closely—both on the training set and extrapolated to the rest of the interval -1 ≤ q ≤ 3.  Clearly, the neural potential trained using 𝒯_3 only captures one well and does not extrapolate correctly to the rest of the domain.

Let V^m(q) denote the neural potential trained on 𝒯_m. We now apply SINDy to the output of each V^m only on its respective training set 𝒯_m, with candidate functions {1, q, q^2, q^3, q^4, q^5, q^6 }.  For reference, the ground truth V(q) can be written as
V(q) = -1 + 2 q + 3 q^2 - 4 q^3 + q^4.
Adjusting λ downward as described above, we find with λ = 0.5 the following algebraic expressions:

    V^1(q)    ≈ -8.138 + 2.0008 q + 3.0009 q^2 - 4.0009 q^3 + 1.0001 q^4 
    V^2(q)    ≈ -6.061 + 2.0032 q + 3.0054 q^2 - 4.0148 q^3 + 0.9748 q^4 
    V^3(q)    ≈ -6.165 + 1.9909 q + 2.9991 q^2 - 3.9886 q^3 + 0.9955 q^4

Noting that the constant terms are irrelevant, we note here that all three models agree closely with the ground truth.  The agreement between V and the algebraic forms of V^1 and V^2 was expected.  We find it somewhat surprising that SINDy, when applied to the output of V^3 on its training set 𝒯^3, yields a quartic polynomial with two wells.



 §.§ Central Force Problem

We consider a central force problem for one particle (d=3) with Hamiltonian

    H(,) = ^2/ 2 + ^-1 + (10 - )^-1.

The norm here is the standard Euclidean norm.  Using explicit Runge-Kutta integration in Mathematica, we generate R=1 trajectory with random initial condition ((0), (0)) chosen uniformly on [-1,1]^6.  Using this trajectory, we compute r(t) = (t) as well as ṙ(t) = dr/dt.  We save the (r(t),ṙ(t)) trajectories at N = 20001 points with a time step of 0.001.  We then search for a reduced-order (d=1) model with Hamiltonian

    H(r,ṙ) = ṙ^2 /2 + V(r),

where V(r) is a neural potential.   We take the neural potential model to have two hidden layers, each with 16 units.  We train for 500000 steps at a learning rate of 10^-3, first using exponential linear unit activations ψ(x) =  x     x ≥ 0 
exp(x) - 1     x < 0.. We then initialize the neural network using the learned weights/biases and retrain using softplus activations ϕ(x) = log(1 + exp(x))—this activation was chosen to enable series expansions of V(r), as described in greater detail below.  Prior to retraining, we also change the network by adding an exponential function to the output layer—we incorporate this function to better model the steep gradients in the potential near r=0 and r=10.  When we retrain, we take 500000 steps at a learning rate of 10^-3.  We carry out the training in two stages because training directly with softplus activations and exponential output failed.



For this problem, classical physics gives us an effective potential

    V_eff(r) = r^-1 + (10-r)^-1 + ℓ^2/(2 r^2).

where ℓ is a conserved quantity determined from the initial condition.  In Figure <ref>, we plot both the trained neural potential V (in red) and the effective potential V_eff(r) (in black).  After adjusting for the constant bias term, we find close agreement.

We then exported the weight and bias matrices to Mathematica, forming the neural potential model

    V(r) = W_3 ϕ( W_2 ϕ(W_1 r + b_1) + b_2) + b_3.

Unlike ψ, the softplus activation ϕ is amenable to series expansion via symbolic computation.   In particular, since we can see that the effective potential V_eff(r) is a rational function, we explored Padé expansions of V(r).  These attempts were unsuccessful in the sense that we did not obtain models of V(r) that are any more interpretable than the compositional form of (<ref>). 

Turning to SINDy, we formed a library of candidate functions

    {1, r^-1, r^-2, r^-3, (10-r)^-1, (10-r)^-2, (10-r)^-3}.

Adjusting λ in the same manner described above, we find that with λ = 0.15, the estimated model is

    V(r) ≈β_0 + β_1 r^-1 + β_2 r^-2 + β_4 (10-r)^-1

Here β_0 ≈ -0.2384, β_1 ≈ 1.005, β_2 ≈ 0.4461, and β_4 ≈ 0.9723.  We see from the form of V_eff(r) given above that β_1 and β_4 are both close to the ground truth values of 1.  Note that for the trajectory on which the system was trained, we have ℓ^2 /2 ≈ 0.4655.  Hence β_2 has an error of less than 4.2%.  This demonstrates a successful application of SINDy to interpret the neural potential as a rational function; this interpretation of V is itself close to V_eff.



 §.§ Charged Particles in Coulomb Potential

We now consider two oppositely charged particles (d=6) subject to the classical Coulomb electrostatic potential.   We take the mass matrix to be M = diag(1, 1/2).  The kinetic energy is T() = ^T M^-1 /2.  If we partition = (_1, _2) where _i is the position of the i-th particle, then the potential is

    V() = -1/4 π1/_1 - _2.

Here we apply the Störmer-Verlet algorithm, a symplectic method, to generate R = 1000 trajectories, each with N = 10001 points recorded at a time step of 0.001.  Each trajectory starts with random initial conditions ((0), (0)) chosen from a standard normal.  For this problem, our goal is to use the data to recover V.  We train two different neural potential models with increasing levels of prior domain knowledge:


  * We first set up the neural network's input layer to compute from  the difference _1 - _2 ∈ℝ^3; the neural potential then transforms this three-dimensional input into a scalar output.  The neural network here has 8 hidden layers, each with 16 units and tanh activations.  Using only 800 of the R = 1000 trajectories, we first train using the first 100 points from each of the 800 trajectories, taking 500000 steps at a learning rate of 0.01. Again restricting ourselves to the 800 training trajectories, we then use the next 100 points, followed by the next 100 points, etc., each time taking 500000 steps at a learning rate of 0.01.  As the training loss was observed to be sufficiently small (≈ 0.006021), we halted training.



In Figure <ref>, we plot both training (left) and test (right) results.  The training results are plotted with the first 1000 points of the 800 trajectories used for training, while the test results are plotted with the first 1000 points of the 200 held out trajectories.  For both plots, we subtracted the maximum computed value of V (on each respective data set).  In each plot, we plot V (on all points  in the training and test sets) versus r = _1 - _2.

Overall, we see reasonable agreement between the neural potential and the ground truth.  Note that the neural network is essentially tasked with discovering that it should compute the inverse of the norm of  _1 - _2.  We suspect that this function of _1 - _2 may be somewhat difficult to represent using a composition of activation functions and linear transformations as in (<ref>).  Despite training for a large number of steps, there is noticeable variation in neural potential values for large r.



We now apply SINDy to V (on the training set) using candidate functions {1, r^-1, r^-2, r^-3}.  Adjusting λ as described above, we find with λ = 0.04, the approximation

    V(r) ≈β_0 + β_1 r^-1

with β_0 ≈ 0.7602 and β_1 ≈ -0.06911.  For comparison, the ground truth coefficient of r^-1 is -(4 π)^-1≈ -0.07958.


  * We then rearchitect the network to include a layer that takes the input  and computes the norm of the difference _1 - _2; the rest of the neural potential is then a scalar function of this scalar input.  Here the neural network has 8 hidden layers, each with 8 units and tanh activations.  We train for 50000 steps with learning rate of 0.05.  Note that here, for training, we use N = 5001 time steps of only 100 trajectories.




In Figure <ref>, we plot both training (left) and test (right) results.  The training results are plotted with the first 5001 points of the 100 trajectories used for training, while the test results are plotted with a completely different set of 100 trajectories, each of length 5001  For both plots, we subtracted the minimum computed value of V (on each respective data set).  In each plot, we compute V on all points  in the training and test sets, and then plot these V values versus r = _1 - _2.

To generate an interpretable version of V (on the training set), we apply SINDy with candidate functions {1, r^-1, r^-2, r^-3}.  Adjusting λ as described above, we find with λ = 0.05, the approximation

    V(r) ≈β_0 + β_1 r^-1

with β_0 ≈ 2.267 and β_1 ≈ -0.07792.  This computed value of β_1 is less than 2.1% away from the ground truth value of -(4 π)^-1; the error for the earlier approximation (<ref>) was just over 13.1%.

Incorporating prior knowledge that the potential should depend only on r dramatically improves the quality of the learned potential.  Essentially, we have eliminated the need for the neural network to learn the norm function.  We outperform the results from Figure <ref> using a less complex network, trained for fewer steps and a larger learning rate.  Comparing with Figure <ref>, we see that Figure <ref> features reduced variation in V(r) for large r, and improved test set results as well.





§ CONCLUSION

We conclude that, for the examples we have explored, our approach does lead to accurate potentials that can themselves be approximated closely by interpretable, closed-form algebraic expressions.
In ongoing/future work, we plan to apply the techniques described here to high-dimensional systems for which reduced-order (i.e., effective) potentials are unknown.  We also seek to extend our method to quantum Hamiltonian systems.  While we have focused here on clean data from known models, we are also interested in learning potentials from noisy time series.  We expect that by adapting the method of <cit.>, we will be able to simultaneously filter the data and estimate an interpretable neural potential.
 







splncs04

10




Artrith2016
Artrith, N., Urban, A.: An implementation of artificial neural-network
  potentials for atomistic materials simulations: Performance for TiO2.
  Computational Materials Science  114,  135–150 (Mar 2016).
  10.1016/j.commatsci.2015.11.047

Behler2016
Behler, J.: Perspective: Machine learning potentials for atomistic simulations.
  The Journal of Chemical Physics  145(17),  170901 (2016).
  10.1063/1.4966192

Behler07
Behler, J., Parrinello, M.: Generalized neural-network representation of
  high-dimensional potential-energy surfaces. Phys. Rev. Lett.  98,
  146401 (Apr 2007). 10.1103/PhysRevLett.98.146401

bhat_nonparametric_2016
Bhat, H.S., Madushani, R.W.M.A.: Nonparametric Adjoint-Based Inference
  for Stochastic Differential Equations. In: 2016 IEEE International
  Conference on Data Science and Advanced Analytics (DSAA). pp.
  798–807 (2016). 10.1109/DSAA.2016.69

Brunton2016
Brunton, S.L., Proctor, J.L., Kutz, J.N.: Discovering governing equations from
  data by sparse identification of nonlinear dynamical systems. Proceedings of
  the National Academy of Sciences  113(15),  3932–3937 (2016)

BhatDale2018
Dale, R., Bhat, H.S.: Equations of mind: Data science for inferring nonlinear
  dynamics of socio-cognitive systems. Cognitive Systems Research  52,
   275–290 (2018)

DaleBhatSindyr
Dale, R., Bhat, H.S.: sindyr: Sparse Identification of Nonlinear Dynamics
  (2018), <https://CRAN.R-project.org/package=sindyr>, r package version
  0.2.1

Duncker2019
Duncker, L., Bohner, G., Boussard, J., Sahani, M.: Learning interpretable
  continuous-time models of latent stochastic dynamical systems. In: Chaudhuri,
  K., Salakhutdinov, R. (eds.) Proceedings of the 36th International Conference
  on Machine Learning. Proceedings of Machine Learning Research, vol. 97, pp.
  1726–1734. PMLR, Long Beach, California, USA (09–15 Jun 2019),
  <http://proceedings.mlr.press/v97/duncker19a.html>

Hansen15
Hansen, K., Biegler, F., Ramakrishnan, R., Pronobis, W., von Lilienfeld, O.A.,
  Muller, K.R., Tkatchenko, A.: Machine learning predictions of molecular
  properties: Accurate many-body potentials and nonlocality in chemical space.
  The Journal of Physical Chemistry Letters  6(12),  2326–2331
  (2015). 10.1021/acs.jpclett.5b00831

Ramakrishnan15
Ramakrishnan, R., Hartmann, M., Tapavicza, E., von Lilienfeld, O.A.: Electronic
  spectra from TDDFT and machine learning in chemical space. The Journal of
  Chemical Physics  143(8),  084111 (2015). 10.1063/1.4928757

BhatICML2019
Raziperchikolaei, R., Bhat, H.S.: A block coordinate descent proximal method for
  simultaneous filtering and parameter estimation. In: Chaudhuri, K.,
  Salakhutdinov, R. (eds.) Proceedings of the 36th International Conference on
  Machine Learning. Proceedings of Machine Learning Research, vol. 97, pp.
  5380–5388. PMLR, Long Beach, California, USA (09–15 Jun 2019),
  <http://proceedings.mlr.press/v97/raziperchikolaei19a.html>

Sahoo2018
Sahoo, S., Lampert, C., Martius, G.: Learning equations for extrapolation and
  control. In: Dy, J., Krause, A. (eds.) Proceedings of the 35th International
  Conference on Machine Learning. Proceedings of Machine Learning Research,
  vol. 80, pp. 4442–4450. PMLR, Stockholmsmässan, Stockholm Sweden (10–15
  Jul 2018), <http://proceedings.mlr.press/v80/sahoo18a.html>

ZhangSchaeffer2018
Zhang, L., Schaeffer, H.: On the Convergence of the SINDy Algorithm.
  arXiv e-prints arXiv:1805.06445 (May 2018)



