On learning Hamiltonian systems from data

Current data science exploration of dynamics often focuses on regression or classification problems, without routinely incorporating additional assumptions about the nature of the system that generated the data. This has started to change recently, with approaches to extract generic dynamical systems by Ref. 32, specifying the variables and possible expressions for the formulas beforehand. In particular, treating the central object to be modeled as the discrete-time flowmap $Φτ(x)=x+∫0τf(x(t))dt$ but learning $Φ$ directly as a black box may result in qualitative differences from the true system.³¹ Using the associated differential equation $dx/dt=f(x)$ instead allows us to exploit established numerical integration schemes to help approximate the flow map and can be accomplished with a neural network by structuring the loss function similarly to such classical numerical integration schemes. In addition to our older work along these lines,^6,25–31 more recent work has revived this approach, focused on treating the layers of a deep neural network as iterates of a dynamical system, where “learning” consists of discovering the right attractors.^3,5,14

In particular, interest has focused on how residual networks¹¹ and highway networks³³ can be interpreted as iterative solvers⁸ or as iterated dynamical systems.⁴ In the latter paper (a NeurIPS 2018 best paper awardee), the authors choose not to unroll the iteration in time explicitly, but instead use continuous-time numerical integration. While the focus was on the dynamics-of-layers concept, time series learning was also performed.

The Koopman operator has also been employed in combination with neural networks to extract conservation laws and special group structures.^12,15 Symmetries in relation to conserved quantities are a well-studied problem in physics.^10,13,19,20 A recent thread of research consists of learning physics models from observation data,¹⁷ including modeling discrete-time data as observations of a continuous-time dynamical system.^21,22

It is informative to study physical systems through their conserved quantities, such as total energy, which can be encoded in a Hamiltonian function.^1,10 The measure-preserving property of Hamiltonian systems has recently been exploited for transport of densities in Markov-chain Monte-Carlo methods¹⁸ and variational autoencoders.^2,24 For our purposes, a natural progression of the thread of computational modeling of physics from observations is to represent the Hamiltonian function directly.

Concurrently with this submission, two papers independently addressing similar issues appeared as preprints.^9,34 In the first of these,⁹ a loss function very similar to parts of our Eq. (12) was used. The second paper focuses on the transformation of densities generated through the Hamiltonian. It is possible to draw some analogies between the second preprint and the old (non-Hamiltonian) work we mentioned above,^6,25–31 as this newer work also used rollout with a time step templated on a classical numerical integration method (here, symplectic Euler and leapfrog). Both papers also use the pendulum as an example, emphasizing conditions where the system can be well approximated by a linear system: a thin annulus around the stable steady state at $(q,p)=(0,0)$ where the trajectories are nearly circular.

The remainder of the paper is structured as follows:

1. We derive data-driven approximations (through two approaches: Gaussian processes and neural networks) of a Hamiltonian function on a given phase space, from time series data. The Hamiltonian functions we consider do not need to be separable as a sum $H(q,p)=T(p)+V(q)$⁠, and in our illustrative example, we always work in the fully nonlinear regime of the pendulum.

2. We build data-driven reconstructions of a phase space from (a) linear and (b) nonlinear, nonsymplectic transformations of the original Hamiltonian phase space. The reconstruction then leads to a symplectomorphic copy of the original Hamiltonian system.

3. We construct a completely data-driven pipeline combining (a) the construction of an appropriate phase space and (b) the approximation of a Hamiltonian function on this new phase space, from nonlinear, high-dimensional observations (e.g., from movies/sequences of movie snapshots).

A Hamiltonian system on Euclidean space $E=R2n$⁠, $n∈N$ is determined through a function $H:E→R$ that defines the equations

where $(⋅):=d/dt$ and $q(t),p(t)∈Rn$ are interpreted as “position” and “momentum” coordinates in the “phase space” $E$⁠. In many mechanical systems, and in all examples we discuss in this paper, the interpretation of the coordinates $q,p$ is reflected in the dynamics through $q˙=p$⁠, i.e., $H(q,p)=12p2+h(q)$ for some function $h:Rn→Rn$⁠. In general, Eqs. (1) and (2) imply that the Hamiltonian is constant along trajectories $(q(t),p(t))$ because

Equations (1) and (2) can be restated as a partial differential equation for $H$ at every $(q,p)∈E$⁠,

where $I∈Rn×n$ is the identity matrix and $ν$ is the vector field on $E$ [the left hand side of (1) and (2)], which only depends on the state $(q,p)$⁠. The symplectic form on the given Euclidean space takes the form of the matrix $ω$⁠.

In Sec. III, we discuss how to approximate the function $H$ from given data points $D={(qi,q˙i,q¨i)}i=1N$⁠. This involves solving the partial differential equation (4) for $H$⁠. Since these equations determine $H$ only up to an additive constant, we assume that we also know the value $H0=H(q0,p0)$ of $H$ at a single point $(q0,p0)$ in the phase space. This is not a major restriction for the approach because $H0$ as well as $(q0,p0)$ can be chosen arbitrarily.

As an example, consider the case $n=1$ and the Hamiltonian

This Hamiltonian forms the basis for the differential equations of the nonlinear pendulum, $q¨=−sin⁡(q)$⁠, or, in the first-order form, $q˙=∂H(q,p)/∂p=p$ and $p˙=−∂H(q,p)/∂q=−sin⁡(q)$⁠. In this section, we numerically solve partial differential equation (PDE) (4) by approximating the solution $H$ using two approaches: Gaussian processes²³ (Sec. III A) and neural networks (Sec. III B).

We model the solution $H$ as a Gaussian process $H^$ with a Gaussian covariance kernel,

where $x$ and $x′$ are points in the phase space, i.e., $x=(q,p)$⁠, $x′=(q′,p′)$⁠, and $ϵ∈R+$ is the kernel bandwidth parameter (we chose $ϵ=2$ in this paper). Given a collection $X$ of $N$ points in the phase space, as well as the function values $H(X)$ at all points in $X$⁠, the conditional expectation of the Gaussian process $H^$ at a new point $y$ is

where we write $k(X,X′)i,j:=k(xi,xj)$ for the kernel matrix evaluated over all $x$ values in the given dataset $X$⁠. In Eq. (7), the dimensions of the symbols are $y∈R2n$⁠, $k(X,X′)∈RN×N$⁠, $k(y,X)∈RN$⁠, and $H(X):=(H(x1),H(x2),…,H(xN))∈RN$⁠. All vectors are column vectors. Estimates of the solution $H$ to the PDE at new points depend on the value of $H$ over the entire dataset. We do not know the values of $H$⁠, but differentiating Eq. (7) allows us to set up a system of equations to estimate $H$ at an arbitrary number $M$ of new points $y$ close to the points $x∈X$⁠, by using the information about the derivatives of $H$ given by the time derivatives $q˙$ in Eq. (1) and (the negative of) $p˙$ in Eq. (2). Let $g(xi)∈R2n$ be the gradient of $H$ at a point $xi$ in the phase space, which we know from data $(q˙,p˙)$⁠. Then,

Together with an arbitrary pinning term at a point $x0$⁠, the list of known derivatives leads to a linear system of $2nN+1$ equations, where we write

for the derivative of the given kernel function $k$ with respect to its first argument, evaluated at a given point $xi$ in the first argument and all points in the dataset $Y$ in the second argument. For each $xi$⁠, we thus have a (column) vector of $M$ derivative evaluations, which is transposed and stacked into a large matrix to form the full linear system

We reiterate: $X$ is a dataset of $N$ points where we know the derivatives of $H$ through $g(xi)=(∂H∂q(xi),∂H∂p(xi))T=(−p˙i,q˙i)T∈R2n$⁠. We evaluate on a fine grid $Y$ of $M$ points (such that $k(Y,Y′)∈RM×M$⁠, $∂∂xk(xi,Y)∈R2n×M$⁠) and have information $g(X)∈R2nN$ on a relatively small set of $N$ points called $X$ [black dots in Fig. 1(b)]. The derivative of the Gaussian process can be stated using the derivative of the kernel $k$ with respect to the first argument. The matrix inverse for $k(Y,Y′)$ is approximated by the inverse of $(k(Y,Y′)+σ2I)$⁠, that is, through Tikhonov regularization with parameter $σ=10−5$⁠, as is the standard for Gaussian process regression. Solving this system of equations for $H(Y)$ yields the approximation for the solution to the PDE. Figure 1(b) shows the result, and Table I lists the mean squared error to the 625 training data points and an independently drawn set of 200 validation points in the same domain, where no derivatives were available. See Ref. 22 for a more detailed discussion of the solution of PDE with Gaussian processes.

Another possibility for learning the form of $H$ using data is to represent the function with an artificial neural network⁷ (ANN). We write

where the activation function $σl$ is nonlinear (except where otherwise indicated, we used $tanh$⁠) for $l=1,…,L$ (if $L≥1$⁠) and the identity for $l=L+1$⁠. The learnable parameters of this ANN are ${(Wl,bl)}l=1,…,L+1$⁠, and we gather all such learnable parameters from the multiple layers that may be used in one experiment into a parameter vector $w$⁠. If there are no hidden layers (⁠$L=0$⁠), then we learn an affine transformation $x1=x0⋅W+b$⁠. This format provides a surrogate function $H^(q,p)=xL+1$⁠, where the input $x0$ is the row vector $[q,p]$⁠. (Treating inputs as row vectors and using right-multiplication for weight matrices is convenient, as a whole batch of $N$ inputs can be presented as an $N$-by-$2$ array.) For all the experiments shown here, this network for computing $H^$ has two hidden layers of width 16.

Similarly, in the case(s) we need to learn additional transformations $θ^$ and $θ^−1$ (see Sec. IV), they are also learned using such networks.

We collect training data by sampling a number of initial conditions in the rectangle $(q,p)∈[−2π,2π]×[−6,6]$⁠, then simulate short trajectories from each to their final $(q,p)$ points. For each of these, we additionally evaluate $(q˙,p˙)$⁠. Shuffling over simulations once per epoch, and dividing this dataset into batches, we then perform batchwise stochastic gradient descent to learn the parameters $w$ using an Adam optimizer on the objective function defined below.

For this paper, all neural networks were constructed and trained using TensorFlow, and the gradients necessary for evaluating the Hamiltonian loss terms in Eq. (12) were computed using TensorFlow’s core automatic differentiation capabilities.

This objective function comprises a scalar function evaluated on each data 4-tuple $d=(q,p,q˙,p˙)$ in the batch and then averaged over the batch. This scalar function is written as

where the dependence on $w$ is through the learned Hamiltonian $H^$⁠, the loss-term weights $ck$ are chosen to emphasize certain properties of the problem thus posed, and the partial derivatives of $∂H^∂p$ and $∂H^∂q$ are computed explicitly through automatic differentiation. Except for $c2$⁠, all $ck$ values are set to either $1$ or $0$ depending on whether the associated loss term is to be included or excluded. Because of the square term in Eq. (5), we set $c2$ arbitrarily to $10$ if nonzero, so the loss is not dominated by $f1$⁠. An alternative might be to set $c1$ to $1/10$⁠.

Since Eqs. (1) and (2) together imply (3), any one of the three terms $f1$⁠, $f2$⁠, and $f4$ could be dropped as redundant; therefore, we can set $c4$ to zero but monitor $H^˙$ as a useful sanity check on the accuracy of the learned solution. In Fig. 1(c), we show the results of this process with our default nonzero values for $ck$⁠.

As an ablation study, we explored the effect (not shown here) of removing the first, second, and fourth terms. By construction, the true $Ht(q,p)$ function is zero for all $(q,p)$⁠. Note that this is only ever achieved to any degree in the central box in the figure, where data were densely sampled. Removing $f4$ made no visible difference in the quality of our $H^t≈0$ approximation, which was expected due to the redundancy in the set of Eqs. (1)–(3). However, removing either $f1$ or $f2$ gives poor results across the figure, despite the apparent redundancy of these terms with $f4$⁠. This might be due to not balancing the contributions of the $p˙$ and $q˙$ terms, for which we attempted to compensate by unequal weighting values $c1$ and $c2$⁠.

We now consider a set of observation functions $θ:E→RM$⁠, $θ=(y1,…,yM)$⁠, with $M≥dim⁡E=2n$⁠, such that $θ$ is a diffeomorphism between the phase space $E$ and its image $θ(E)$⁠. In this setting, the notion of a symplectomorphism is important.¹ In general, a symplectomorphism is a diffeomorphism that leaves the symplectic structure on a manifold invariant. In our setting, a symplectomorphism of $E=Q×P$ maps to a deformed space $E^=Q^×P^$ where the system dynamics in the new variables $q^∈Q^$⁠, $p^∈P^$ is again Hamiltonian and conjugate to the original Hamiltonian dynamics. Not every diffeomorphism is a symplectomorphism, and we do not assume that $θ$ is a symplectomorphism. A constant scaling of the coordinates is also not possible to distinguish from a scaling in the Hamiltonian function itself, so the recovered system will be a symplectomorphic copy of the original, scaled by an arbitrary constant.

In the setting of this section, we do not assume access to $E$⁠, $H$⁠, or the explicit form of $θ$⁠. Only a collection of points $θi$ and time derivatives $ddtθi$ in the image $θ(E)$ is available. We describe an approach to approximate a new map $θ^−1:θ(E)→E^$ into a symplectomorphic copy of $E$ through an autoencoder,⁷ such that the transformed system in $E^$ is conjugate to the original Hamiltonian system in $E$⁠. Upon convergence, and if we had access to $θ$⁠, the map $θ^−1°θ≡S:E→E^$ would approximate a symplectomorphism, and $θ−1°θ^≡S−1$ would be its inverse. During the estimation of $θ^−1$⁠, we simultaneously approximate the new Hamiltonian function $H^:E^→R$⁠. Figure 2(a) visualizes the general approach, where only the information $(x,y)i$ and $ddt(x,y)i$ is available to the procedure, while $θ^−1$ and a Hamiltonian $H^$ on $E^$ are constructed numerically.

An interesting and important feature common to the next three examples is that one cannot expect to systematically recover the original $(q,p)$ values from the given observation data $(x,y)$⁠. Only symplectic transformations $(q^,p^)$ can be recovered, which are enough to define a Hamiltonian. Once the coordinates $(q^,p^)$ are fixed, the Hamiltonian function in these coordinates is unique up to an additive constant.

Table I contains training and validation loss of the networks for all experiments. Additionally, we trained the network from Sec. IV D with only 331 images and did observe a significantly higher validation loss (not shown), consistent with overfitting.

The following loss function is used to train an autoencoder component together with a Hamiltonian function approximation network component:

where the dependence on $w$ is through the learned Hamiltonian $H^$ and the learned transformations $θ^$ and $θ^−1$⁠, and the time derivatives in the space $E^$ are computed as $q^˙=x˙∂q^∂x+y˙∂q^∂y$ and $p^˙=x˙∂p^∂x+y˙∂p^∂y$⁠, using the Jacobian $Dθ^−1=∂q^∂x∂q^∂y∂p^∂x∂p^∂y$ of the transformation $θ^−1$ (computed pointwise with automatic differentiation). When we learn a (especially nonlinear) transformation $θ^−1$⁠, in addition to the Hamiltonian $H^(q^,p^)$⁠, including the $f4$ term in the composite loss can have a detrimental effect on the learned transformation. There exists an easily-encountered naïve local minimum in which $θ^−1$ maps all of the sampled values from $θ(E)∋(x,y)$ to a single point in $E^$⁠, and the Hamiltonian learned is merely the constant function at the pinning value, $H^(q^,p^)=H^0$⁠. In this state (or an approximation of this state), all of $∂H∂q^$⁠, $∂H∂p^$ and the elements of $Dθ^−1$ are zero, so the loss (with terms $f1$⁠, $f2$⁠, $f3$⁠, and optionally $f4$⁠) is zero (resp., small). A related failure is that in which the input in $θ(E)$ is collapsed by $θ^−1$ to a line or curve in $E^$⁠.

To alleviate both of these problems, we added a new loss component $f6$⁠. That is, we require that the learned transformation does not collapse the input. It is sufficient for the corresponding weighting factor $c6$ to be a very small nonzero value (e.g., $10−6$⁠). The addition of $f6$ to our loss helps us to avoid falling early in training into the unrecoverable local minimum described above and also helps keep the scale of the transformed variables $q^$ and $p^$ macroscopic.

We generate data from a rectangular region $x∈[−1,1]$⁠, $y∈[−1,1]$ and then transform the region linearly with $θ−1(x,y)=A−1[x,y]T=[q,p]T$⁠. The matrix $A−1$ is the inverse of $A=R⋅Λ$⁠; a scaling is followed by a rotation where $Λ=λ100λ2$⁠, $λ1=1$⁠, $λ2=64$⁠, $R=cos⁡ρ−sin⁡ρsin⁡ρcos⁡ρ$⁠, and $ρ=5°$⁠. Our observation data $θ(E)$ are thus given by $[x,y]T=A⋅[q,p]T$⁠. Using the true Hamiltonian $H(q,p)=p2/2+(1−cos⁡q)$⁠, we additionally compute true values for $dq/dt$ and $dp/dt$ and then use $A$ to propagate these to $x$ and $y$ via $dxdt=∂x∂qdqdt+∂y∂pdpdt$ and similar for $y$⁠, where the partial derivatives are computed analytically (here, just the elements of $A$ itself).

Our network is then presented with observation data $x,y$ and its corresponding time-derivatives. Its task is to learn $A^$ and $A^−1$⁠, which convert to and from variables $q^,p^$ (symplectomorphic to the original $q,p$⁠) and a Hamiltonian $H^$ in this new space. When evaluating the loss, the time derivatives of $q^$ and $p^$ are likewise computed via automatic differentiation using the chain rule through the learned transformation, e.g., as $dq^dt=∂q^∂xdxdt+∂q^∂ydydt$⁠. Note also that $[q^,p^]T=A^−1⋅A⋅[q,p]T$⁠, so if the original space $E$ could be found, $A^$ would satisfy $A^⋅A−1=I$⁠. This cannot be expected given only the data in $θ(E)$⁠; we can only be sure that $A^⋅A−1$ approximates a symplectomorphism of the original $E$⁠.

We could learn $θ^−1$ from a general class of nonlinear functions, as a small $tanh$ neural network, but here, we simply learn $A^$ and $A^−1$ as linear transformations (that is, we have a linear “neural network,” where $L=0$ in Eq. (11), and $b1$ is fixed as $0$⁠). As we include the reconstruction error of this autoencoder in our loss function, $A^−1$ is constrained to be the inverse of $A^$ to a precision no worse than the $f1$ and $f2$ terms in Eq. (14), after all three are scaled by their corresponding $ck$ values. In fact, for the linear case, initially the autoencoder’s contribution to the loss is significantly lower than the Hamiltonian components (see Fig. 3), but, as training proceeds and the $f1$ and $f2$ terms are improved (at the initial expense of raising the autoencoder loss), larger reductions in loss are possible by optimizing $H^$ rather than $θ^$⁠, so $f5$ is decreased as quickly as (the larger of) $f1$ or $f2$⁠. That is, the autoencoder portion of the loss falls quickly to the level where it no longer contributes to the total loss given its weighting in the loss sum.

We find that the learned symplectomorphism $S(q,p)=A^−1⋅A⋅[q,p]T$⁠, depicted in its $q$ portion in Fig. 3, preserves $q$ unmixed with $p$ in one or the other of its two discovered coordinates. This is because both (a) $(q,p)↦(p,−q)$ as well as (b) $(q,p)↦(q,p+f(q))$ for any smooth function $f$ are symplectomorphisms. They are special because $H(q,p)=H^(q^(q,p),p^(q,p))$⁠; i.e., they even preserve the Hamiltonian formulation. For the map (a), the transformation of the Hamiltonian can be seen from the following derivation:

Here, the first equality of (15) and (16) follows from the map, and the last equality of (15) and (16) follows from the requirement that $q^,p^$ follow Hamiltonian dynamics with respect to the new Hamiltonian $H^$⁠. Equations (18) and (20) then show that the new Hamiltonian is the same as the old one (modulo an additive constant) when considered as a map on the old coordinates.

In addition to the linear $θ$ of Sec. IV B, we show comparable results for a nonlinear transformation $θ$ and learned $θ^$⁠. Specifically, we transform the data through $(x,y)=θ(q,p)$ where

the inverse of which is given by $q=(x−y2)20$ and $p=(y−x2+2xy2−y4)10$⁠. We use the analytical Jacobian of this $θ$ to compute the necessary $x˙$ and $y˙$ for input to our network.

We proceed as before, except that we no longer restrict the form of the learned $θ^$ and $θ^−1$ to linear transformations, but instead allow small multilayer perceptrons of a form similar to that used for $H^$⁠.

The resulting induced symplectomorphism is again one that appears to preserve an approximately monotonically increasing or decreasing $q$ in either $q^$ or $p^$⁠. This can be seen in Fig. 4.

As a further demonstration of the method, we use a graphical rendering of the moving pendulum example from before as the transformation $θ$ from the intrinsic state $(q,p)$ to an image $x$ as our high-dimensional observable. We use a symplectic semi-implicit Euler’s method,

to generate $q(ti),p(ti)$ trajectories for various initial conditions and then a simple graphical renderer to display these as images [see Fig. 5(a)]. When rendering our video frames, we drag a tail of decaying images behind the moving pendulum head so that information about both position $q$ and velocity $p$ is present in each rendered frame. This is done by iterating over each $q(ti),p(ti)$ trajectory and, for each $ti$⁠, (1) multiplying the entire current image by a cooling factor of $∼0.96$⁠, (2) adding a constant heating amount to the image in a circle of fixed radius centered around the current $cos⁡(q),sin⁡(q)$ point, and (3) clipping the image per-pixel to lie within $[0,1]$⁠. Samples of the resulting images are visible in Fig. 5(a).

Though we do not use $p(ti)$ directly in this procedure, its value is observed in the length of the resulting tail dragged behind the moving pendulum head. We create trajectories long enough that the effect of the initial formation of the tail (during which its length is not necessarily a good indicator of $p$⁠) is not visible anymore and then use only the final two observations from these trajectories.

In order to make the approach agnostic to the data, we do not want to assume that the space $E^$ is periodic, so instead, we use a four-dimensional phase space with elements $z^=[q^1,q^2,p^1,p^2]=[q^,p^]$ and consider the splitting into $(q^1,q^2)$ and $(p^1,p^2)$ during training. In the space of input images, the manifold does not fill up four-dimensional space but a cylinder, which is mapped to the four-dimensional encoding layer by the autoencoder.

In addition, to simplify the learning problem, we learn $θ^−1$ as the combination of a projection onto the first 20 principal components of the training dataset followed by a dense autoencoder, reserving learning $θ^−1$ as an end-to-end convolutional autoencoder for future work. The encoding provides $z^$ and, as before, we learn $θ^−1$ in tandem with $H^$⁠, where now the conditions of Eq. (14) are upgraded to vector equivalents to accommodate $z^$⁠.

In Sec. IV A, we added a loss term proportional to the reciprocal of the determinant of the transformation’s Jacobian in order to avoid transformations that collapsed the phase space. Here, this was not such an issue—of course, some collapsing of the high-dimensional representation is obviously required. Instead, a common mode of failure turned out to be learning constant $H^$ functions, which automatically satisfy the Hamiltonian requirements (a constant is naturally a conserved quantity). To avoid this, we considered several possible ways to promote a nonflat $H^$ function, ultimately settling on (a) adding a term that encouraged the standard deviation of $H^$ values to be nonzero and (b) minimizing not just the mean squared error in our $f1$ and $f2$ terms, but also the max squared error, to avoid trivial or bilevel $H^(q,p)$ functions. The spherical Gaussian prior used in training variational autoencoders³⁴ also avoids, as a by-product, learning constant functions.

The result, shown in Fig. 5(b), was a pulled-back $H^(q,p)$ function that at least in broad strokes resembles the truth and that satisfies $dH^/dt≈0$ (typically about $10−2$⁠).

We described an approach to approximate Hamiltonian systems from observation data. It is a completely data-driven pipeline to (a) construct an appropriate phase space and (b) approximate a Hamiltonian function on the new phase space, from nonlinear, possibly high-dimensional observations (here, movies).

When only transformations of the original Hamiltonian phase space can be observed, it is only possible to recover a symplectic copy of the original phase space, with additional freedom in a constant scaling of the coordinates and an additive constant to the Hamiltonian function. If no additional information about the original space is available, this is a fundamental limitation and not particular to our approach. It is not necessary that there is an “original phase space” at all, and so, the resulting symplectic phase space is by no means unique. In fact, Darboux's theorem implies the opposite.^35,36 There is no local obstruction to a symplectic structure, so all symplectic manifolds of dimension 2n are locally symplectomorphic to the standard flat space $R2n$⁠. The choice of a single such space has to rely on other factors, such as, possibly, interpretability by humans or simplicity of the equations.

The approach may be extended to time-dependent Hamiltonian functions. This would allow us to cope with certain dissipative systems.¹⁶ An even broader extension may allow transformations to arbitrary normal forms as the “target vector field” and thus would not be constrained to Hamiltonian systems. In the general case, it will become important to explore whether the transformation we approximate remains bounded over our data or whether it starts showing signs of approaching a singularity, suggesting that the problem may not be solvable.

This work was funded by the US Army Research Office (ARO) through a Multidisciplinary University Research Initiative (MURI) and by the Defense Advanced Research Projects Agency (DARPA) through their Physics of Artificial Intelligence (PAI) program, and through agreement HR001118C0100 (I.M.).