Learning and Interpreting Potentials for Classical Hamiltonian Systems

TL;DR

This paper presents a method to learn and interpret potential energy functions from trajectories of classical Hamiltonian systems, combining neural networks with equation discovery to produce interpretable models.

Contribution

The paper introduces a novel approach that learns neural potentials and extracts closed-form algebraic expressions, enabling interpretable modeling of classical Hamiltonian systems.

Findings

Accurately learns neural potentials matching ground truth.

Successfully extracts interpretable algebraic expressions.

Demonstrates effectiveness on multiple classical systems.

Abstract

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.

Equations

$$H(\mathbf{p},\mathbf{q})=T(\mathbf{p})+V(\mathbf{q}).$$

$$\dot{\mathbf{q}}$$

$$\dot{\mathbf{p}}$$

$$L(\mathbf{\theta })=\frac{1}{RN}\sum _{j=1}^R\sum _{i=0}^{N-1}{\Vert \frac{{\mathbf{p}}_{i+1}^j-{\mathbf{p}}_i^j}{h}+{\nabla }_{\mathbf{q}}V({\mathbf{q}}_i^j;\mathbf{\theta })\Vert }^2.$$

$$\mathbf{V}=\Xi \mathbf{\beta }+ϵ$$

$$\hat{V}(\mathbf{x})=\sum _{j=1}^J{\beta }_j{\xi }^j(\mathbf{x})+ϵ.$$

$$H(q,p)=\frac{p^2}{2}+\frac{q^2}{2}.$$

$$\hat{V}(q)\approx {\beta }_0+{\beta }_2{q}^2$$

$$V(q)=x^2(x-2{)}^2-(x-1{)}^2.$$

$${\hat{V}}^1(q)$$

$${\hat{V}}^2(q)$$

$${\hat{V}}^3(q)$$

$$H(\mathbf{q},\mathbf{p})=\frac{\Vert \mathbf{p}{\Vert }^2}{2}+\Vert \mathbf{q}{\Vert }^{-1}+(10-\Vert \mathbf{q}\Vert {)}^{-1}.$$

$$H(r,\dot{r})=\frac{{\dot{r}}^2}{2}+\tilde{V}(r),$$

$$V_{\text{eff}}(r)=r^{-1}+(10-r{)}^{-1}+{\ell }^2/(2r^2).$$

$$\hat{V}(r)=W_3\varphi (W_2\varphi (W_{1}r+b_1)+b_2)+b_3.$$

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

$$\hat{V}(r)\approx {\beta }_0+{\beta }_1{r}^{-1}+{\beta }_2{r}^{-2}+{\beta }_4(10-r{)}^{-1}$$

$$V(\mathbf{q})=-\frac{1}{4\pi }\frac{1}{\Vert {\mathbf{q}}_1-{\mathbf{q}}_2\Vert }.$$

$$\hat{V}(r)\approx {\beta }_0+{\beta }_1{r}^{-1}$$

Full text

1 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\geq 1$, let $\mathbf{q}(t)\in\mathbb{R}^{d}$ and $\mathbf{p}(t)\in\mathbb{R}^{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

[TABLE]

In this paper, we consider the problem of learning or identifying the potential energy $V$ from data $(\mathbf{q}(t),\mathbf{p}(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 $(\widetilde{\mathbf{q}}(t),\widetilde{\mathbf{p}}(t))\in\mathbb{R}^{2\widetilde{d}}$. The direct approach is to integrate numerically the $6m$-dimensional system of differential equations for the full Hamiltonian (1) and then use the full solution $(\mathbf{q}(t),\mathbf{p}(t))$ to compute $(\widetilde{\mathbf{q}}(t),\widetilde{\mathbf{p}}(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 $(\widetilde{\mathbf{q}}(t),\widetilde{\mathbf{p}}(t))$ itself satisfies a Hamiltonian system, we can search for a potential $\widetilde{V}(\widetilde{\mathbf{q}})$ that yields an accurate, reduced-order model for $(\widetilde{\mathbf{q}}(t),\widetilde{\mathbf{p}}(t))$. If $\widetilde{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 $(\widetilde{\mathbf{q}}(t),\widetilde{\mathbf{p}}(t))$ directly, with computational savings that depend on $d/\widetilde{d}$.

There is a rapidly growing literature on machine learning of potential energies in computational/physical chemistry, e.g., [3, 9, 10, 2, 1]. 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., [5, 4, 6, 12, 8], 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.

2 Approach

Assume $T(\mathbf{p})=\sum_{i=1}^{d}M_{ii}^{-1}p_{i}^{2}$ where $M$ is a diagonal mass matrix. Then, from (1), we can write Hamilton’s equations:

[TABLE]

Let the training data consist of a set of $R$ trajectories; the $j$-th such trajectory is $\{\mathbf{q}^{j}_{i},\mathbf{p}^{j}_{i}\}_{i=0}^{N}$. Here $(\mathbf{q}^{j}_{i},\mathbf{p}^{j}_{i})$ denotes a measurement of $(\mathbf{q}(t),\mathbf{p}(t))$ at time $t=ih$ 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 (2a). We now posit a model for $V$ that depends on a set of parameters $\boldsymbol{\theta}$. For instance, if we model $V$ using a neural network, $\boldsymbol{\theta}$ stands for the collection of all network weights and biases. Then we use (2b) to form an empirical risk loss

[TABLE]

Let $\tau=Nh$ denote the final time in our grid. Note that (3) approximates $E\left[(1/\tau)\int_{t=0}^{t=\tau}\|\dot{\mathbf{P}}(t)+\nabla_{\mathbf{Q}}V(\mathbf{Q}(t);\boldsymbol{\theta})\|^{2}\mathrm{dt}\right]$, the expected mean-squared error of a random trajectory $(\mathbf{Q}(t),\mathbf{P}(t))$ assumed to satisfy (2a).

We model $V$ using a dense, feedforward neural network with $L\geq 2$ layers. Because we train with multiple trajectories, the input layer takes data in the form of two tensors—one for $\mathbf{q}$ and one for $\mathbf{p}$—with dimensions $N\times R\times d$. The network then transposes and flattens the data to be of dimension $NR\times 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 $\nu$ units where $1\leq\nu\leq d$. As these architectural details differ by example, we give them below.

Note that the loss (3) involves the gradient of $V$ with respect to the input $\mathbf{q}$. We use automatic differentiation to compute this gradient. More specifically, in our IPython/Jupyter notebooks (linked below), we use the batch_jacobian method in TensorFlow. This is easy to implement, fast, and accurate up to machine precision.

The trained network gives us a neural potential $\widehat{V}:\mathbb{R}^{d}\to\mathbb{R}$. To interpret $\widehat{V}$, we apply the SINDy method [5]. We now offer a capsule summary of this technique. Suppose we have a grid $\{\mathbf{x}^{k}\}_{k=1}^{K}$ of points in $\mathbb{R}^{d}$. We use the notation $\mathbf{x}^{k}=(x_{1}^{k},\ldots,x_{d}^{k})$. We evaluate $\widehat{V}$ on the grid, resulting in a vector of values that we denote by $\mathbf{V}$. We also evaluate on the grid a library of $J$ candidate functions $\xi^{j}:\mathbb{R}^{d}\to\mathbb{R}$ for $1\leq j\leq J$; each such evaluation results in a vector $\Xi^{j}$ that we take to be the $j$-th column of a matrix $\Xi$. In $d=1$, examples of candidate functions are $\{1,x,x^{2},\ldots\}$ or $\{1,x^{-1},x^{-2},\ldots\}$. In $d=2$, an example is $\{1,x_{1},x_{2},x_{1}^{2},x_{1}x_{2},x_{2}^{2},\ldots\}$. Each candidate function is simply a scalar-valued function on $\mathbb{R}^{d}$.

Equipped with the $K\times 1$ vector $\mathbf{V}$ and the $K\times J$ matrix $\Xi$, we solve the regression problem

[TABLE]

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

As shown recently [13], this algorithm converges in a finite number of steps to an approximate minimizer of $\|\mathbf{V}-\Xi\boldsymbol{\beta}\|^{2}+\lambda^{2}\|\boldsymbol{\beta}\|_{0}$. Here $\|\boldsymbol{\beta}\|_{0}$ denotes the number of nonzero entries of $\boldsymbol{\beta}$. Hence, increasing the parameter $\lambda$ leads to a more sparse set of coefficients $\boldsymbol{\beta}$. Once we have fit the regression model in this way, we obtain an interpretable model of $\widehat{V}$, specifically:

[TABLE]

If $\boldsymbol{\beta}$ is highly sparse, most of the coefficients $\beta_{j}$ will be zero. Suppose that the candidate functions $\xi^{j}$ are well-known functions such as positive or negative powers of the coordinates $x_{i}$ of the input $\mathbf{x}$. 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 $\epsilon$ here captures the error in this sparse approximation of $\widehat{V}$. In general, one chooses $\lambda$ to balance the sparsity of $\boldsymbol{\beta}$ with the quality of the approximation $\|\epsilon\|$.

3 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 sindyr package [7] in R to interpret the neural potential. In what follows, the mass matrix $M$ in (2) 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.

3.1 Simple Harmonic Oscillator

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

[TABLE]

Exact trajectories consists of circles centered at the origin in $(q,p)$ space. For training data, we use $R=10$ such circles; for $1\leq i\leq 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 $\widehat{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 1, we plot both the trained neural potential $\widehat{V}$ (in red) and the true potential $V(q)=q^{2}/2$ (in black). When plotting $\widehat{V}$, we have subtracted a constant bias (the minimum obtained value of $\widehat{V}$) so that the curve reaches a minimum value of zero. Note that this constant bias is completely unimportant for physics, as only $\nabla V$ appears in Hamilton’s equations (2) and the loss (3). However, the constant bias does lead us to include an intercept in the regression model (4), 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 $\widehat{V}(q)$ with candidate functions $\{1,q,q^{2},q^{3}\}$. In the following test, and in fact throughout this paper, we start with $\lambda=1$ and tune $\lambda$ downward until the error $\|\epsilon\|$—between the neural network potential $\widehat{V}(q)$ and the SINDy-computed approximation—drops below $10^{-10}$. We find that with $\lambda=0.04$, the estimated system is

[TABLE]

with $\beta_{0}\approx-49.18$ and $\beta_{2}\approx 0.4978$. We see that $\widehat{V}$ closely tracks the true potential $V(q)=q^{2}/2$ up to the constant bias term, which can be ignored.

3.2 Double Well

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

[TABLE]

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 $\mathcal{T}_{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 $\mathcal{T}_{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 $\mathcal{T}_{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 $\mathcal{T}_{m}$, we train for $50000$ steps at a learning rate of $0.01$.

We seek to understand how the choice of training set $\mathcal{T}_{m}$ affects the ability of the neural potential $\widehat{V}$ to track the true potential (8). We plot and discuss the results in Figures 2 and 3. Overall, the neural potentials trained using $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$ match $V(q)$ closely—both on the training set and extrapolated to the rest of the interval $-1\leq q\leq 3$. Clearly, the neural potential trained using $\mathcal{T}_{3}$ only captures one well and does not extrapolate correctly to the rest of the domain.

Let $\widehat{V}^{m}(q)$ denote the neural potential trained on $\mathcal{T}_{m}$. We now apply SINDy to the output of each $\widehat{V}^{m}$ only on its respective training set $\mathcal{T}_{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+2q+3q^{2}-4q^{3}+q^{4}$. Adjusting $\lambda$ downward as described above, we find with $\lambda=0.5$ the following algebraic expressions:

[TABLE]

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 $\widehat{V}^{1}$ and $\widehat{V}^{2}$ was expected. We find it somewhat surprising that SINDy, when applied to the output of $\widehat{V}^{3}$ on its training set $\mathcal{T}^{3}$, yields a quartic polynomial with two wells.

3.3 Central Force Problem

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

[TABLE]

The norm here is the standard Euclidean norm. Using explicit Runge-Kutta integration in Mathematica, we generate $R=1$ trajectory with random initial condition $(\mathbf{q}(0),\mathbf{p}(0))$ chosen uniformly on $[-1,1]^{6}$. Using this trajectory, we compute $r(t)=\|\mathbf{q}(t)\|$ as well as $\dot{r}(t)=dr/dt$. We save the $(r(t),\dot{r}(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

[TABLE]

where $\widetilde{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 $\psi(x)=\begin{cases}x&x\geq 0\\ \exp(x)-1&x<0.\end{cases}$. We then initialize the neural network using the learned weights/biases and retrain using softplus activations $\phi(x)=\log(1+\exp(x))$—this activation was chosen to enable series expansions of $\widehat{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

[TABLE]

where $\ell$ is a conserved quantity determined from the initial condition. In Figure 4, we plot both the trained neural potential $\widehat{V}$ (in red) and the effective potential $V_{\text{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

[TABLE]

Unlike $\psi$, the softplus activation $\phi$ is amenable to series expansion via symbolic computation. In particular, since we can see that the effective potential $V_{\text{eff}}(r)$ is a rational function, we explored Padé expansions of $\widehat{V}(r)$. These attempts were unsuccessful in the sense that we did not obtain models of $\widehat{V}(r)$ that are any more interpretable than the compositional form of (12).

Turning to SINDy, we formed a library of candidate functions

[TABLE]

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

[TABLE]

Here $\beta_{0}\approx-0.2384$, $\beta_{1}\approx 1.005$, $\beta_{2}\approx 0.4461$, and $\beta_{4}\approx 0.9723$. We see from the form of $V_{\text{eff}}(r)$ given above that $\beta_{1}$ and $\beta_{4}$ are both close to the ground truth values of $1$. Note that for the trajectory on which the system was trained, we have $\ell^{2}/2\approx 0.4655$. Hence $\beta_{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 $\widehat{V}$ is itself close to $V_{\text{eff}}$.

3.4 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=\operatorname{diag}(1,1/2)$. The kinetic energy is $T(\mathbf{p})=\mathbf{p}^{T}M^{-1}\mathbf{p}/2$. If we partition $\mathbf{q}=(\mathbf{q}_{1},\mathbf{q}_{2})$ where $\mathbf{q}_{i}$ is the position of the $i$-th particle, then the potential is

[TABLE]

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 $(\mathbf{q}(0),\mathbf{p}(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 $\mathbf{q}$ the difference $\mathbf{q}_{1}-\mathbf{q}_{2}\in\mathbb{R}^{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 ($\approx 0.006021$), we halted training.

In Figure 5, 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 $\widehat{V}$ (on each respective data set). In each plot, we plot $\widehat{V}$ (on all points $\mathbf{q}$ in the training and test sets) versus $r=\|\mathbf{q}_{1}-\mathbf{q}_{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 $\mathbf{q}_{1}-\mathbf{q}_{2}$. We suspect that this function of $\mathbf{q}_{1}-\mathbf{q}_{2}$ may be somewhat difficult to represent using a composition of activation functions and linear transformations as in (12). Despite training for a large number of steps, there is noticeable variation in neural potential values for large $r$.

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

[TABLE]

with $\beta_{0}\approx 0.7602$ and $\beta_{1}\approx-0.06911$. For comparison, the ground truth coefficient of $r^{-1}$ is $-(4\pi)^{-1}\approx-0.07958$. 2. 2.

We then rearchitect the network to include a layer that takes the input $\mathbf{q}$ and computes the norm of the difference $\|\mathbf{q}_{1}-\mathbf{q}_{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 6, 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 $\widehat{V}$ (on each respective data set). In each plot, we compute $\widehat{V}$ on all points $\mathbf{q}$ in the training and test sets, and then plot these $\widehat{V}$ values versus $r=\|\mathbf{q}_{1}-\mathbf{q}_{2}\|$.

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

[TABLE]

with $\beta_{0}\approx 2.267$ and $\beta_{1}\approx-0.07792$. This computed value of $\beta_{1}$ is less than $2.1\%$ away from the ground truth value of $-(4\pi)^{-1}$; the error for the earlier approximation (15) 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 5 using a less complex network, trained for fewer steps and a larger learning rate. Comparing with Figure 5, we see that Figure 6 features reduced variation in $\widehat{V}(r)$ for large $r$, and improved test set results as well.

4 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 [11], we will be able to simultaneously filter the data and estimate an interpretable neural potential.