Discovering conservation laws using optimal transport and manifold learning

Conservation laws are key theoretical and practical tools for understanding, characterizing, and modeling nonlinear dynamical systems. However, for many complex systems, the corresponding conserved quantities are difficult to identify, making it hard to analyze their dynamics and build stable predictive models. Current approaches for discovering conservation laws often depend on detailed dynamical information or rely on black box parametric deep learning methods. We instead reformulate this task as a manifold learning problem and propose a non-parametric approach for discovering conserved quantities. We test this new approach on a variety of physical systems and demonstrate that our method is able to both identify the number of conserved quantities and extract their values. Using tools from optimal transport theory and manifold learning, our proposed method provides a direct geometric approach to identifying conservation laws that is both robust and interpretable without requiring an explicit model of the system nor accurate time information.


Introduction
Conservation laws are powerful constraints on the dynamics of many physical systems in nature, and the corresponding conserved quantities are essential features for characterizing the behavior of these systems.Through Noether's theorem, conservation laws are closely tied with the symmetries of a physical system and play a key role in our understanding of physics.Conservation laws also help stabilize and enhance the performance of predictive models for complex nonlinear dynamics, e.g.symplectic integrators for Hamiltonian systems [1] and pressure projection for incompressible fluid flow [2].In fact, for chaotic dynamical systems, conserved quantities are often the only features of the system state that can be reliably known far into the future.Discovering conservation laws helps us characterize the long term behavior of complex dynamical systems and understand the underlying physics.
While the conservation laws of many physical systems are well-known and often derived from known symmetries, there are still many instances where it is difficult to even determine the number of conservation laws, let alone explicitly extract the conserved quantities.As a historical example, consider the Korteweg-De Vries (KdV) equation modeling shallow water waves.The KdV equation, despite its apparent complexity, has infinitely many conserved quantities [3] and is, in fact, fully solvable via an inverse scattering transform [4]a discovery made after significant theoretical and computational effort.Developing better general methods for identifying conserved quantities will allow us to improve our understanding of new or understudied physical systems and build more efficient and stable predictive models.
In real-world applications, an accurate model for the underlying physical system is often unavailable, forcing us to identify conservation laws using only sample trajectories of the system dynamics.One broad approach is to use modern data-driven methods based on the Koopman operator formulation of dynamical systems, which lifts the dynamics into an infinite dimensional operator space [5].In the Koopman formalism, conserved quantities are just one type of Koopman eigenfunction with eigenvalue zero.Thus, one approach is to first apply a system identification method, such as dynamic mode decomposition [6,7], sparse identification with a library of basis functions [8], or even deep learning-based approaches [9][10][11], to model the system dynamics and then set up and solve the Koopman eigenvalue problem.Alternatively, previous work has also proposed directly setting up the eigenvalue problem by estimating time derivatives from data and then fitting the conserved quantities using a library of possible terms [12] or a neural network [13].These methods can work quite well but require that the measured trajectories have sufficiently low noise and high time resolution in order to accurately estimate time derivatives.
Constructing a model for a dynamical system provides much more information than just the conservation laws.In fact, even estimating time derivatives is usually not necessary if we are only interested in identifying conserved quantities.In this work, we will instead focus on an alternative approach that does not require an explicit model or detailed time information but rather takes advantage of the geometric constraints imposed by conservation laws.Specifically, the presence of conservation laws restricts each trajectory in phase space to lie solely on a lower dimensional isosurface of the conserved quantities.The dimensionality of these isosurfaces can provide information about the number of conserved quantities or constraints [14].Furthermore, since each isosurface corresponds to a particular set of conserved quantities, the variations in shape of the isosurfaces directly correspond to variations in the conserved quantities.In other words, we can identify and extract conserved quantities by examining the varying shapes of the isosurfaces sampled by the trajectories.
In contrast with recent work using black box deep learning methods to fit conserved quantities that are consistent with the sampled isosurfaces [15,16], we propose and demonstrate a nonparametric manifold learning approach (Fig. 1) that directly characterizes the variations in the sampled isosurfaces, producing an embedding of the space of conserved quantities.Our method first uses the Wasserstein metric from optimal transport [17] to compute distances in shape space between pairs of sampled isosurfaces and then extracts a low dimensional embedding for the manifold of isosurfaces using diffusion maps [18,19].Each point in this embedding corresponds to a distinct isosurface and therefore to a distinct set of conserved quantities, i.e. the embedding explicitly parameterizes the space of varying conserved quantities.Related methods have been recently suggested for characterizing molecular conformations using the 1-Wasserstein distance together with diffusion maps [20], performing system identification by comparing invariant measures using the 2-Wasserstein distance [21], and reconstructing normal forms using diffusion maps [22].Recent theoretical work has also formalized the idea of using alternative non-Euclidean norms, like the Wasserstein distance, in spectral embedding methods such as diffusion maps [23].
We provide an analytic analysis of our approach for a simple harmonic oscillator system and numerically test our method on several physical systems: the single and double pendulum, planar gravitational dynamics, the KdV equation for shallow water waves, and a nonlinear reactiondiffusion equation that generates an oscillating Turing pattern.We also demonstrate the robustness of our approach to noise in the measured trajectories, missing information in the form of a partially observed phase space, as well as approximate conservation laws (additional experiments in Appendix E and Appendix C).In our comparison tests (Appendix F), our approach outperforms prior deep learning-based direct fitting methods, all while being an order of magnitude faster.We also provide an easy-to-use codebase, which parallelizes across multiple GPUs, to make an efficient implementation of our method as accessible as possible. 1 Proposed Manifold Learning Approach Our proposed approach uses manifold learning to identify and embed the manifold of phase space isosurfaces sampled by the trajectories of a dynamical system.In particular, we compute a diffusion map over a set of trajectories, each of which samples a particular phase space isosurface (Fig. 1a).The pairwise distances between these trajectories are given by the 2-Wasserstein distance (Fig. 1b), providing the metric structure necessary for applying diffusion maps (Fig. 1c).
The manifold embedding extracted by the diffusion map corresponds directly to the space of conserved quantities (Fig. 1d).Note that this type of analysis does not require knowledge of the equations of motion (Eq. 1) and makes no direct reference to time.

Dynamical Systems
Consider a dynamical system with states x ∈ M that live in a d-dimensional phase space M and evolve in time according to a system of first order ODEs dx dt = F(x) with n conserved quantities G 1 (x), . . ., G n (x).

Conserved Quantities and Phase Space Isosurfaces
A conserved quantity G i (x) is a function of the system state x that does not change over time, i.e.
but may vary across different initial conditions.As a result, along a particular trajectory x(t), the n conserved quantities form a set of constraints which depend on the values of the conserved quantities c = {c 1 , c 2 , . . ., c n }.This set of constraint equations restricts the trajectory to lie in a phase space isosurface X c ⊆ M with dimension d − n.In fact, if any point of a trajectory lies on the isosurface X c , then all other points from the trajectory will lie on the same isosurface.By studying the variations in shape of these isosurfaces, we are able to directly characterize the space of conserved quantities.In particular, consider the manifold C formed by the isosurfaces X c in shape space. 2 This manifold C is parameterized by the conserved quantities c.Therefore, by analyzing C using manifold learning, we can extract the conservation laws of the underlying dynamical system.

Ergodicity and Physical Measures
To uniquely identify the isosurface associated with each trajectory, we must make several additional assumptions that will allow us to treat the set of points making up each trajectory as samples from an ergodic invariant measure on the corresponding isosurface.Specifically, we assume that, for each trajectory x(t) with conserved quantities c, the dynamical system (Eq. 1) admits a physical measure [24] that is ergodic on the isosurface X c and is defined by where δ x(t) is the Dirac measure centered on x(t).This ensures that trajectories with the same conserved quantities will sample the same distribution on the same isosurface, allowing us to use the distribution sampled by each trajectory as a proxy for the corresponding isosurface.For more details about this assumption, see Appendix G.4.
In practice, the sampled distribution may be lower dimensional than the corresponding isosurface if some of the conserved quantities do not vary in the dataset and instead correspond to fixed constraints, or if the dynamical system is dissipative.In the former case, this does not affect our  1 Proposed non-parametric method for discovering conservation laws illustrated using a simple pendulum example.(a) First, we collect and normalize the trajectory data from the dynamical system.Two example trajectories are highlighted in red and blue.(b) Then, we use the Wasserstein metric from optimal transport to compute the distance between each pair of trajectories and construct a distance matrix.For the two example trajectories, the optimal transport plan is shown as lines connecting pairs of points.The constructed distance matrix is plotted with color representing the computed Wasserstein distance between each pair of trajectories.The computed distance between the two example trajectories is marked (black dots) on the distance matrix plot.(c) An embedding of the shape space manifold C is extracted from the distance matrix using diffusion maps.The embedding plot is colored by the conserved energy of the pendulum E. The points corresponding to the two example trajectories are marked in red and blue.(d) Finally, a heuristic score (Appendix A) is used to select relevant components.In this case, only component 1 is relevant, corresponding to a single conserved quantity-the energy E. Again, the embedding plot is colored by E, and the two example trajectories are marked in red and blue.
ability to uniquely identify a distribution with an isosurface and its corresponding set of conserved quantities, meaning that we are able to apply this approach even if the provided phase space is much larger than the intrinsic phase space of the dynamical system.In the latter case, the dissipative nature of the system may cause information about conservation laws relevant during the transient portion of the dynamics to be lost, but we are still able to use our approach to identify conserved quantities relevant for the long term behavior of the system.

Wasserstein Metric
To analyze the isosurface shape space manifold C-i.e. the manifold of conserved quantitiesusing manifold learning methods, we need to place some structure on the points X c ∈ C. Having associated each isosurface X c with a corresponding distribution defined by an ergodic physical measure µ c , we choose to lift the Euclidean metric on the phase space into the space of distributions using the 2-Wasserstein metric from optimal transport where the cost function c(x, y) = ∥x − y∥ 2 is the squared Euclidean distance, and π ∈ Π(µ c , µ c ′ ) is a valid transport map between µ c and µ c ′ [17].
For discrete samples, the 2-Wasserstein distance between two sets of sample points {x 1 , x 2 , . . ., x S } and {y 1 , y 2 , . . ., y S } is defined as where the cost matrix C ij = ∥x i − y j ∥ 2 , and the transport matrix T is subject to the constraints To efficiently compute an entropy regularized form of this optimization problem, we use the Sinkhorn algorithm [25] and estimate the Wasserstein distance as a debiased Sinkhorn divergence [26].
One important subtlety in this construction is the choice of the ground metric for optimal transport.As previously mentioned, we use a Euclidean metric on the phase space, which implicitly imposes a choice of units to make the phase space dimensionless.In fact, there is no canonical choice for the ground metric, and different choices result in different Wasserstein metrics on the shape space.While, in theory, information about all conserved quantities will be embedded in the resulting distance matrix regardless of the choice of metric, the metric ultimately determines how easy it is to access this information.For example, when multiple conserved quantities are present, the relative effect of each conserved quantity on the computed Wasserstein distances will determine how prominent each conserved quantity is and how easily it is identified using manifold learning.To partially mitigate this issue and improve consistency, we normalize each component of our data to have a maximum absolute value of 1 before computing the pairwise Wasserstein distances.
Finally, using the Wasserstein distance provides our approach with a tremendous amount of robustness (Appendix E), but also makes it susceptible to certain kinds of sampling inhomogeneity.See Appendix G.3 for a more detailed discussion of this trade off.

Diffusion Maps
Using the structure provided by the Wasserstein metric, we then use diffusion maps to generate an embedding for C. The diffusion map manifold learning method uses a spectral embedding algorithm applied to an affinity matrix to construct a low dimensional embedding of the data manifold [18,19].Using the pairwise Wasserstein distances W 2 (µ i , µ j ) computed from discrete samples provided by the trajectory data (Eq.6), we first construct a kernel matrix using a Gaussian kernel and then scale it to form an affinity matrix for our spectral embedding where D i = k K ik , and α is a hyperparameter.The spectral embedding algorithm then takes this affinity matrix and constructs a normalized graph Laplacian where I is the identity matrix.The eigenvectors v i corresponding to the smallest eigenvalues λ i ≥ 0 (excluding λ 0 = 0) of the Laplacian then provide an approximate low dimensional embedding of the manifold of conserved quantities C. In our experiments, we set α = 1 so that the Laplacian computed by the spectral embedding algorithm approximates the Laplace-Beltrami operator [18].
To estimate the dimensionality of C and to choose which eigenvectors v i to include in our embedding, we use a heuristic score that combines a measure of relevance, given by a length scale computed from the Laplacian eigenvalues, with a previously suggested measure of "unpredictability" for minimizing redundancy [27] (alternative approaches also exist [28,29]).To construct our embedding, we only include the Laplacian eigenvectors with score above a chosen cutoff value and discard the rest as either noise or redundant embedding components.To determine the cutoff, we perform a sweep of the cutoff value looking for robust ranges and find that a cutoff of 0.6 works well across all of our experiments, which consist of a wide variety of datasets and dynamical systems.See Appendix A for more details.

Analytic Result for the Simple Harmonic Oscillator
In the case of a simple harmonic oscillator (SHO) without measurement noise and in the infinite sample limit, we are able to explicitly derive an analytic result for our proposed procedure.We first compute the pairwise distances provided by the Wasserstein metric and then derive the embedding produced by a diffusion map, which corresponds to the conserved energy of the SHO.

Wasserstein Metric: Constructing the Isosurface Shape Space
Consider a SHO with Hamiltonian given in terms of position q and momentum p.
The SHO energy isosurfaces E = H(q, p) form concentric ellipses in a 2D phase space.Choosing units such that m = 1 and ω = 1, we obtain concentric circles with uniformly distributed samples (assuming a uniform sampling in time).The 2-Wasserstein distance between a pair of uniformly distributed circular isosurfaces is simply given by the difference in radii |r 1 −r 2 |.This is because, due to the rotational symmetry of the two distributions, the optimal transport plan for an isotropic cost function is to simply move each point on isosurface 1 radially outward (or inward) to the point on isosurface 2 with the same angle θ.This result does not meaningfully change with a different choice of units, which is equivalent to rescaling the phase space coordinates q, p.If we rescale q, p by factors k q , k p , our cost function simply becomes where we label points on the isosurfaces by their angle θ on the original circular isosurfaces.The SHO optimal transport plan Π takes θ on isosurface 1 to the point with the same angle θ on isosurface 2, and Π for the SHO is invariant to coordinate rescaling (Appendix I).Therefore, the total transport cost is so the 2-Wasserstein distance is i.e. the same result modulo a constant factor.While this is not a general result, we find that our approach is often fairly robust to such changes, including the extreme case of scaling some phase space coordinates all the way down to zero resulting in a partially observed phase space (Appendix E).

Diffusion Maps: Extracting the Conserved Energy
Once we have pairwise distances in the isosurface shape space, we can use diffusion maps to study the resulting manifold of isosurface shapes.With sufficient samples, the operator constructed by the diffusion map should converge to the Laplace-Beltrami operator on the manifold.For the SHO, the isosurface shape space is isomorphic to R + with each circular isosurface mapped to its radius.
If we sample trajectories with radii r ∈ (0, √ 2E 0 ) for some maximum energy E 0 , then the manifold is a real line segment, and the resulting Laplacian operator (with open boundary conditions) has eigenvalues λ n = π 2 n 2 /2E 0 and corresponding eigenvectors v n (r) = cos( √ λ n r).Therefore, the first eigenvector or embedding component successfully encodes the conserved energy and is, in fact, a monotonic function of the energy.

Numerical Experiments
To demonstrate and empirically test our method for discovering conservation laws, we generate datasets from a wide range of dynamical systems, each consisting of randomly sampled trajectories with different initial conditions and the corresponding conserved quantities.Note that we use the dimensionless form of each dynamical system.All of the code necessary for reproducing our results is available at https://github.com/peterparity/conservation-laws-manifold-learning.

Simple Harmonic Oscillator
We first numerically test our analytic result for the SHO and obtain good agreement (Fig. 2) using both the default scaling k q = k p = 1 (Figs.2ad) as well as the position only scaling k q = 1, k p = 0 (Figs.2e-h), which effectively reduces the dimension of the phase space.A linear fit of the first embedding component from the diffusion map with the analytically predicted component  (Eq.15) achieves a correlation coefficient of R 2 = 0.9995 for the default scaling and R 2 = 0.9961 for the position only scaling.We also verify that the heuristic score (Appendix A) accurately determines that there is only one relevant embedding  component (Figs.2c,g), which corresponds to the conserved energy.

Simple Pendulum
To demonstrate our method on a simple nonlinear dynamical system, we analyze a simple pendulum that has a 2D phase space consisting of the angle θ and angular momentum ω (Fig. 3a).The equations of motion are This system has a single scalar conserved quantity corresponding to the total energy of the pendulum, so the trajectories form 1D orbits in phase space (Fig. 3b).
Our method is able to correctly determine that there is only a single conserved quantity (Fig. 3c) corresponding to the energy of the pendulum (Fig. 3d).The single extracted embedding component is monotonically related to the energy with Spearman's rank correlation coefficient ρ = 0.9997.We are also able to achieve similar results (ρ = 0.9978) with a high level of Gaussian noise (standard deviation σ = 0.5) added to the raw trajectory data (Figs.3e-h), showing that our approach is quite robust to measurement noise.

Planar Gravitational Dynamics
To test our method on a system with multiple conserved quantities, we simulate the gravitational system of a planet orbiting a star with much greater mass (Fig. 4a).We fix the orbits to all lie in a 2D plane, giving us an effectively 4D phase space.The resulting equations of motion are This system has one scalar and two vector conserved quantities which, in our 4D phase space, reduces to three scalar conserved quantities: the total energy E (or equivalently, the semi-major axis a = −1/2E), the angular momentum L = |L|, and the orbital orientation angle ϕ, which is the angle of the LRL vector A relative to the x-axis.As a result, the trajectories also form 1D orbits in the phase space (Fig. 4b).
Our approach accurately identifies the three conserved quantities (Fig. 4c), and the extracted embedding corresponds most directly to the geometric features of the orbits (Figs.4d-f).The first two components embed the semi-major axis vector a = (a cos ϕ, a sin ϕ) with magnitude given by the semi-major axis a = −1/2E, which is related to the energy E, and orientation given by the orientation angle ϕ of the elliptical orbit (Figs.4d,e).The third relevant component (component 6) embeds the angular momentum L (Fig. 4f).See Appendix A.1 for details on choosing a cutoff to identify the relevant components.A linear fit of the identified relevant embedding components with a cos ϕ (a sin ϕ) has R 2 = 0.987 (R 2 = 0.986) and rank correlation ρ = 0.994 (ρ = 0.992).A similar linear fit with L has R 2 = 0.927 and ρ = 0.970.This example demonstrates that, for a system with multiple conserved quantities, the ground metric for optimal transport controls the relative scale of each conserved quantity in the extracted embedding.In this case, the geometry of the shape space C is dominated by changes in the semi-major axis a and orientation angle ϕ, whereas changes in the angular momentum L, which controls the eccentricity of the orbit, play a more minor role and thus appear in a later embedding component with a lower score (Fig. 4c).

Double Pendulum
To test our approach on a non-integrable system with higher dimensional isosurfaces, we study the classic double pendulum system (Fig. 5a) with unit masses and unit length pendulum arms.This system has a 4D phase space, consisting of the angles θ 1 , θ 2 and the angular velocities ω 1 , ω 2 of the two pendulums (Fig. 5b), and only has a single scalar conserved quantity ) corresponding to the total energy.However, the double pendulum system has both chaotic and Colored by E + + E non-chaotic phases.In particular, at high energies, the system is chaotic and only conserves the total energy, while at low energies, the system behaves more like two coupled harmonic oscillators with two independent (approximately) conserved energies corresponding to the two modes of the coupled oscillator system.Therefore, we expect to see two distinct phases in our extracted embedding: one with a single conserved quantity E at high energy and another with two approximately conserved quantities E ± at low energy, which approximately sum to At first glance, it appears as though our method has only identified a single relevant component corresponding to the conserved total energy E (Figs. 5c,e) with rank correlation ρ = 0.996.However, if we restrict ourselves to low energy trajectories with first embedding component v 1 < −1, we find that there is a region of the shape space that is two-dimensional, corresponding to the two independently conserved energies E ± of the low energy non-chaotic phase where the double pendulum behaves like a coupled oscillator system with two distinct modes.For the low energy trajectories, a linear fit of the now two relevant components with E + (E − ) has rank correlation ρ = 0.919 (ρ = 0.937).If we restrict ourselves to even lower energy trajectories with v 1 < −2, a similar linear fit for E + (E − ) has rank correlation ρ = 0.990 (ρ = 0.989).
This analysis of the double pendulum shows that our method can still provide significant insight into complex dynamical systems with multiple phases involving varying numbers of conserved quantities.This manifests itself as manifolds of different dimensions in shape space that are stitched together at phase transitions, presenting a significant challenge for most manifold learning methods.In this example, this difficulty is reflected in the performance of the heuristic score (Fig. 5c,d), which has trouble telling whether the embedding is one or two dimensional precisely because it is a combination of both a one and two dimensional manifold.The embedding, on the other hand, remains very informative despite the sudden change in dimensionality and allows us to identify interesting features of the system, such as nonlinear periodic orbits (see Appendix D).The effectiveness of diffusion maps when handling these complex situations has been previously observed in parameter reduction applications [30] and is worth studying in more detail in the future.

Oscillating Turing Patterns
Next, we consider an oscillating Turing pattern system that is both dissipative and has a much higher dimensional phase space than our previous examples.In particular, we study the Barrio-Varea-Aragón-Maini (BVAM) model [31,32] with D = 0.08, C = −1.5, and H = 3, following Aragón et al. [32] who showed that this set of parameters results in a spatial Turing pattern that also exhibits chaotic oscillating temporal dynamics, on a periodic domain with size 8.The phase space of the BVAM system consists of two functions u(x) and v(x) 3 which we discretize on a mesh of size 50, giving us an effective phase space dimension of 100.Because this system is dissipative, we will focus on characterizing the long term behavior of the dynamics, i.e. the oscillating Turing pattern, which appears to have a conserved spatial phase η for our chosen set of parameters corresponding to the spatial position of the Turing pattern.In the language of dynamical systems, η parameterizes a continuous set of attractors for this dissipative system.
Our method successfully identifies the spatial phase η but embeds the angle as a circle in a 2D embedding space (Fig. 6)-a result of the periodic topology of η.While this shows that the number of relevant components determined by our heuristic score may not always match the true manifold dimensionality, such cases are often easily identified by examining the components directly (Fig. 6c The heuristic score (with cutoff 0.6) identifies two relevant components, but on further examination, (c) we see that there is just a single conserved angle, corresponding to the spatial phase η of the Turing pattern, that needs to be embedded in two dimensions due to its topology.
an intrinsic dimensionality estimator [33].A linear fit of the two relevant components with cos η (sin η) has R 2 = 0.9991 (R 2 = 0.9997) and ρ = 0.9993 (ρ = 0.9992).This example both tests our method on a high dimensional phase space and demonstrates how our approach can be applied to dissipative systems to study long term behavior.

Korteweg-De Vries Equation
For many spatiotemporal dynamical systems, the conservation laws are local in nature.Locality can significantly simplify the analysis of the conserved quantities and suggests a way to restrict the type of conserved quantities identified by our method.Specifically, we can adapt our approach to focus on local conserved quantities by replacing the raw states (Fig. 7a) by a distribution of local features (Fig. 7b), removing the explicit spatial label and providing a fully translation invariant representation of the state.Then, instead of using the Euclidean metric in the original phase space, we use the energy distance [34,35] between the distributions of local features as the ground metric for optimal transport.To demonstrate this method for identifying local conserved quantities, we consider the Korteweg-De Vries (KdV) equation The KdV equation is fully integrable [4] and has infinitely many conserved quantities [3], the most robust of which are the most local conserved quantities expressible in terms of low order spatial derivatives.To focus on these robust local conserved quantities, we use finite differences (i.e.u(x), ∆u(x) = u(x + ∆x) − u(x), ∆ 2 u(x), . ..) as our local features, allowing us to restrict the spatial derivative order of the identified conserved quantities.In this experiment, we only take u(x) and ∆u(x), meaning that the identified local conserved quantities will only contain up to first order spatial derivatives.For the KdV equation, there are three such local conserved quantities: where c 1 and c 2 are often identified as "momentum" and "energy", respectively [3].These local conserved quantities also have direct analogues in generalized KdV-type equations, hinting at their robustness [36].
Our method successfully identifies three relevant components (Fig. 7c) corresponding to (d-f) the three local conserved quantities (Eq.24).Linear fits of these components to c 1 , c 2 , and c 3 have rank correlations ρ = 0.995, 0.994, and 0.985, respectively.This result shows how our approach can be adapted to incorporate known structure, such as locality and translation symmetry, in applications to complex high dimensional dynamical systems.

Discussion
We have proposed a non-parametric manifold learning method for discovering conservation laws, tested our method on a wide variety of dynamical systems-including complex chaotic systems with multiple phases and high dimensional spatiotemporal dynamics-and also shown how to adapt our approach to incorporate additional structure such as locality and translation symmetry.While our experiments use dynamical systems with known conserved quantities in order to validate our approach, our method does not require any a priori information about the conserved quantities.Our method also does not assume or construct an explicit model for the system nor require accurate time information like previous approaches [12,13], only relying on the ergodicity of the dynamics modulo the conservation laws (Sec.2.1.2).As a result, our approach is also quite robust to measurement noise and can often deal with missing information such as a partially observed phase space (Figs.2e-h, Figs.3e-h, Appendix E).
Compared with recently proposed deep learning-based methods [15,16], our approach is substantially more interpretable since it relies on explicit geometric constructions and well-studied manifold learning methods that directly determine the geometry of the shape space C and, therefore, the identified conserved quantities.This is reflected in our ability to explicitly derive the expected result for the simple harmonic oscillator (Sec.3), as well as in the identified conserved quantities in many of our experiments.For example, the embedding of the semi-major axis vector in the planar gravitational dynamics experiment (Sec.4.3) stems directly from the elliptical geometry of the orbits and their orientation in phase space, which is captured by the Euclidean ground metric and lifted into shape space by optimal transport.Our method also correctly captures the subtleties of the double pendulum system (Sec.4.4) by providing an embedding that shows both a 1D manifold at high energies and a 2D manifold at low energies-a difficult prospect for deep learning approaches that try to explicitly fit conserved quantities.In addition, we empirically find that our method outperforms existing direct fitting approaches [15,16].See Appendix F for a comparison benchmark using our planar gravitational dynamics dataset.
Our manifold learning approach to identifying conserved quantities provides a new way to analyze data from complex dynamical systems and uncover useful conservation laws that will ultimately improve our understanding of these systems as well as aid in developing predictive models that accurately capture long term behavior.While our method does not provide explicit symbolic expressions for the conserved quantities (which may not exist in many cases), we do obtain a full set of independent conserved quantities, meaning that any other conserved quantity will be a function of the discovered ones.Our method also serves as a strong non-parametric baseline for future methods that aim to discover conservation laws from data.Finally, we believe that similar combinations of optimal transport and manifold learning have the potential to be applied to a wide variety of other problems that also rely on geometrically characterizing families of distributions, and we hope to investigate such applications in the near future.
Acknowledgments.We would like to acknowledge useful discussions with Justin Solomon, Ziming Liu Author contributions.All three authors contributed to the conception, design, and development of the proposed method.P.Y.L. implemented, tested, and refined the method and also generated the datasets.The manuscript was written by P.Y.L. with support from R.D. and M.S.

Appendix A Heuristic Score for a Minimal Diffusion Maps Embedding
Traditionally, diffusion maps [19] and Laplacian eigenmaps [18] leave the embedding dimension n as a hyperparameter and simply use the eigenvectors corresponding to the n smallest eigenvalues to construct the embedding.In practice, the embedding dimension n is often chosen for convenience (e.g. in visualization applications) or by examining the eigenvalues λ i and looking for a sharp increase in the magnitude of the eigenvalues that would separate the signal from the noise.Because identifying the number of conservation laws is an important step in our approach, we refine this heuristic by directly computing an approximate length scale where ϵ is the scale factor from the Gaussian kernel (Eq.8) used to construct the Laplacian matrix L. We derive this length scale by considering the normalized kernel I − L to be an approximation of the heat kernel exp(ϵ∆), implying that the length scales l i associated with the Laplace-Beltrami operator ∆ are given by exp(ϵ∆ We then divide by l 1 to obtain the relative length scale which can be used as a heuristic measure of relevance-components with a small relative length scale are more likely to be noise.Compared with directly using the eigenvalues λ i , we find this heuristic to be less sensitive to the choice of ϵ in the kernel.If we examine these three components, we find that they together embed a second order angular mode of components 1 and 2 (Figs.4d,e).In particular, the embedding is shaped like the surface of a cone with the height (or radial distance) roughly corresponding to the semi-major axis a and the angle around the cone corresponding to ϕ mod π, a second order mode of the orientation angle ϕ.
In addition to noise, there is the common problem of redundant embedding components that stem from the structure of the Laplacian operator: higher order modes of previous eigenvectors often appear before more informative eigenvectors corresponding to new manifold directions.This problem is clearly illustrated in the planar gravitational dynamics experiment (Sec.and relevant conserved quantity (Figs.A1d-f).To address this issue, the key observation is that, while all components of the diffusion map are linearly independent, redundant components are still predictable (via a nonlinear function) from previous components.Therefore, we require a measure of "unpredictability" that allows us to identify redundancies.We choose the heuristic m i proposed by Pfau and Burgess [27] that uses a nearest neighbor estimator (using 5 nearest neighbors) to determine whether a new embedding component is too predictable and therefore redundant.Alternative methods for dealing with these redundant components, also called repeated eigendirections or higher harmonics, have been proposed that use local linear regression to detect redundant components [28] or adapt the diffusion kernel to be anisotropic for chosen components [29].
Our final heuristic score (Fig. A1c) is the product of the relative length scale l i /l 1 (Fig. A1a) and the unpredictability measure m i (Fig. A1b).We find this simple combined score performs well for identifying relevant embedding components by removing both noise components as well as redundant components.

A.1 Choosing a Score Cutoff
To use the heuristic score to identify the number of conserved quantities and construct a minimal embedding, we require a score cutoff to separate relevant components that we keep in our embedding from irrelevant components that we discard.
To choose this cutoff, we sweep cutoff values in the interval [0, 1], compute the embedding size (i.e. the number of relevant components) based on the chosen cutoff, and then examine the result to identify a robust value for the cutoff (Fig. A2).Specifically, we look for wide plateaus in the embedding size that indicate robustness to the value of the cutoff and find that a cutoff of 0.6 works well in all of our experiments.In practice, we would treat a cutoff of 0.6 as a good starting point but recommend analyzing a range of cutoff values to find a robust choice of cutoff, as illustrated in Fig. A2.

Appendix B Additional Method Details
All of the code necessary for generating our datasets, applying our method, and reproducing our results is available at https://github.com/peterparity/conservation-laws-manifold-learning.

B.1.1 Entropy Regularized Optimal Transport
We compute an approximate 2-Wasserstein distance using the Sinkhorn algorithm [25], which solves an entropy regularized relaxation of the optimal transport problem ) where the cost matrix C ij = ∥x i − y j ∥ 2 and the entropy h(T ) = − i,j T ij log T ij .W 2 reduces to the exact 2-Wasserstein distance W 2 (Eq.6) as γ → 0. For γ > 0, the entropy regularization introduces a smoothing bias that manifests as a nonzero "self-distance" W 2 ({x i }, {x j }) > 0. This can be corrected by instead using the Sinkhorn divergence [26,35] (B6) as our estimate for the 2-Wasserstein distance, which explicitly subtracts off this self-distance.In our experiments, we use a convergence threshold of ε = 0.01 and a decaying entropy regularization parameter γ that starts at 10.0 and decays by a factor of 0.995 at each step until it reaches a target of 0.1.
Note that the Sinkhorn algorithm is a general approach for solving entropy-regularized optimal transport and is equally good at approximating a 1-Wasserstein distance (or Earth mover's distance).We experimented with using 1-Wasserstein vs. 2-Wasserstein distances and did not find a significant difference in the performance of our method.

B.1.2 Time Complexity
With S samples per trajectory, the Sinkhorn algorithm solves the entropy regularized optimal transport problem in O(S 2 log S/ε 2 ) time for an ε-accurate solution [37,38] using O(S) space (without explicit storing the cost matrix C [39]).Therefore, the time complexity of computing approximate Wasserstein distances for all pairs of trajectories is O(N 2 S 2 log S/ε 2 ) for a dataset containing N total trajectories.This computation is currently the performance bottleneck of our approach (Appendix G.1) but is easily parallelized over multiple GPUs using the OTT-JAX library [39].

B.2.1 Choice of Manifold Learning Method
Unlike many standard manifold learning applications, our input to the manifold learning method is not a set of points in Euclidean space but rather a pairwise Wasserstein distance matrix.Many popular methods, such as local linear embedding [40] and local tangent space alignment [41], explicitly require the data to be embedded in a Euclidean  (b, c) Similarly to the diffusion map embedding, Isomap components 1 and 2 embed the semi-major axis vector a with magnitude a = −1/2E related to the energy and orientation given by the angle ϕ.These components have high rank correlation ρ = 0.994 (ρ = 0.992) with a cos ϕ (a sin ϕ).(d) Component 3 roughly corresponds to the angular momentum L but is noisier than the embedding identified by diffusion maps (Fig. 4f) and has a much lower rank correlation ρ = 0.723 with L.
space and so are not applicable for our problem.In addition to diffusion maps, we also experimented with Isomap [42], which can also take a distance matrix as input.However, we found Isomap embeddings to be noisier and less reliable than diffusion map embeddings (e.g.see Fig. B3).Diffusion maps [19] or Laplacian eigenmaps [18] also provide an effective way to estimate manifold dimensionality and choose relevant embedding components (Appendix A).

B.2.2 Gaussian Kernel Width
The primary hyperparameter in our diffusion maps algorithm is the width parameter ϵ of the Gaussian kernel (Eq.8).Generally speaking, we should choose ϵ to be large enough to avoid noise induced by sparse sampling but small enough for the true heat kernel to be well approximated by the Gaussian kernel [18].We choose ϵ such that the corresponding standard deviation σ = ϵ/2 of the Gaussian kernel is equal to the maximum distance to the kth nearest neighbor.
In practice, especially when the relevant embedding components are well separated from the noise in terms of length scale (Eq.A1), we find the diffusion map to be fairly insensitive to the choice of k, which we generally set to k = 20 nearest neighbors.The only exception is for the planar gravitational dynamics dataset, where we use k = 200.This is because the angular momentum-the least prominent of the three conserved quantities-has a slightly poorer reconstruction for k = 20 (ρ = 0.910) and gets pushed back to component 10 of the embedding.While it is still clearly identifiable from noise, it requires a lower heuristic score cutoff of around ∼ 0.4 and is easier to miss.

B.2.3 Noise Robustness
To improve the noise robustness of our diffusion map, we follow Karoui and Wu [43] and replace the diagonal of the affinity matrix M (Eq.9) with zeros, i.e.
before constructing the Laplacian matrix L.
Because this induces an overall shift in the eigenvalues of the Laplacian that interacts poorly with our length scale heuristic (Eq.A1), we correct for this by subtracting off the normalized mean shift from the Laplacian matrix L to obtain the corrected Laplacian which we use to generate our embeddings.

B.2.4 Time Complexity
For a fixed number of embedding dimensions and a dense kernel matrix K with N × N entries (N being the number of trajectories), diffusion maps have time complexity O(N 2 ) and space complexity O(N 2 ).In our experiments, computing the diffusion map is very fast with run times under one second for every dataset we tested.

B.2.5 Out-of-Sample Embedding
One complication of manifold learning methods like diffusion maps is that they do not provide an explicit way to embed new out-of-sample data.A naive approach would be to rerun the diffusion map algorithm on a combined dataset consisting of the original data used to create the embedding and the new data.However, this does not retain the original embedding and, in some cases, may be computationally prohibitive.A popular approach that does retain the original embedding is the Nyström method [44], which embeds a new point using its pairwise distances with the original data and scales as O(N ).For even faster embedding, landmark diffusion maps offer a significant speed-up by choosing a small subset of M ≪ N landmark points to use during Nyström out-ofsample embedding [45].This also has the added benefit of a reduced memory footprint, since only the M landmark points need to be retained for embedding.

Appendix C Langevin Harmonic Oscillator
To demonstrate our method on a simple example of a dynamical system with approximately conserved quantities at short time scales but no true conserved quantities, we consider an underdamped harmonic oscillator that is weakly coupled to a heat bath.The resulting dynamics are governed by the Langevin equations where the damping γ = 10 −2 .ξ(t) is a Gaussian random process (i.e.Brownian motion) with zero mean and ⟨ξ(t)ξ(t ′ )⟩ = 2γτ δ(t − t ′ ).We set the temperature of the heat bath to be τ = 2 × 10 −2 .At short times t ≪ γ −1 and t ≪ (γτ ) −1/2 , the energy E = (q 2 + p 2 )/2 is approximately conserved.At long times t ≫ γ −1 , all trajectories will sample a stationary distribution with variance ⟨q 2 ⟩ 0 = τ and no conserved quantities.
Fig. C4 Identifying conserved quantities for the Langevin harmonic oscillator over varying sampling times T .(a) Over a short time scale T = 2×10 1 , the timeaveraged energy ⟨E⟩ still clearly distinguishes the different trajectories, so our approach identifies a single (approximately) conserved quantity that corresponds to ⟨E⟩.(b) Component 1 of the embedding is highly correlated with ⟨E⟩ and still matches well with the theoretical result for the simple harmonic oscillator (Eq.15).(c,d) Over a slightly longer time scale T = 2 × 10 2 , we see a very similar result with a noisier fit between component 1 and ⟨E⟩.(e,f) For T = 2 × 10 3 , we start to see more ambiguity, with the identified component becoming significantly less prominent and an even noisier fit.(g,h) Over a long time scale T = 2×10 4 , the system has settled into a stationary distribution and no longer has any even approximately conserved quantities.
We generate four datasets for this system, each with 200 trajectories and 200 samples over different periods of time, i.e. t ∈ [0, T ] with sampling time T ∈ {2 × 10 1 , 2 × 10 2 , 2 × 10 3 , 2 × 10 4 }.This allows us to study how changing the sampling time T -the time scale over which we identify approximately conserved quantities-changes the embedding produced by our approach.Over short time scales, the time-averaged energy ⟨E⟩ is still a distinguishing feature of the trajectories and acts as an approximately conserved quantity (Figs.C4ad).Over longer time scales, the energy becomes less and less relevant until, finally, the system reaches a stationary distribution (Figs.C4e-h).

C.1 Separation of Time Scales and Dependence on γ and τ
As previously noted, the energy of the Langevin harmonic oscillator is only conserved over time scales T much less than the time scales associated with dissipation γ −1 and random forcing (γτ ) −1/2 .However, this approximate conservation law is only meaningful if the energy is roughly conserved over a time scale T much greater than the period of oscillation, which is of order one in our units.
The rank correlation ρ (Fig. C5a) shows the alignment of the first embedding component with the time-averaged energy ⟨E⟩, while the prominence s 1 −max({s i } i>1 ) (Fig. C5b) of the heuristic score shows how well distinguished the first component is from the remaining noisy components.We find that for small γ ≤ 0.01 and τ ≤ 0.02, our approach successfully identifies the approximately conserved energy at intermediate time scales T = 20 or T = 200 (as shown by the high rank correlation and high score prominence).For longer time scales or larger γ and τ , the approximate conservation law no longer holds and so our method identifies no conserved quantities (as shown by the low score prominence).

Appendix D Nonlinear Periodic Orbit of the Double Pendulum
In addition to the chaotic and linear non-chaotic phases, the double pendulum can also exhibit other kinds of complex behavior, including highly nonlinear periodic orbits.In our extracted embedding (Fig. D6a), we see an example of such a nonlinear periodic orbit (Figs.D6b,c).The placement of this periodic orbit in the embedding also meaningfully connects it with the low energy inphase mode from the linear coupled oscillator regime (Fig. 5g), i.e. this periodic orbit can be thought of as a nonlinear high energy extension of the low energy in-phase mode.

Appendix E Robustness to Noise & Partial Observations
To further demonstrate the robustness of our approach, we show several additional experiments on the simple pendulum, planar gravitational dynamics, and double pendulum datasets.For the simple pendulum, our method still performs well when using only angle θ measurements, i.e. a partially observed phase space (Fig. E7a).In fact, even if we add Gaussian noise (standard deviation σ = 0.5) to the raw trajectory data in addition to using a partially observed phase space, we still obtain a similar result (Fig. E7b).Similarly, for planar gravitational dynamics with only position r data or with added Gaussian noise (σ = 0.5), our method is still able to identify the three conserved   E7 Additional experiments illustrating the robustness of our approach.(a) For the simple pendulum system, even when provided only angle θ measurement data (without angular velocity ω), our method is able to identify a single relevant component corresponding to the energy the pendulum (ρ = 0.998).(b) If we then also add σ = 0.5 Gaussian noise, we can still achieve a similar result (ρ = 0.996).For planar gravitational dynamics, our method also performs well given (c) only position r data or (d) with σ = 0.5 Gaussian noise, correctly identifying the three conserved quantities.
quantities (Figs.E7c,d).For the double pendulum, we again see that we retain the same detail in the extracted embedding using only position data (Fig. E8).The corresponding rank correlations with the ground truth conserved quantities are given in Table E1.
The robustness of our method to both noise and partial observations is largely a consequence of using the Wasserstein distance as our metric for comparing trajectories.In our problem formulation, we consider the isosurfaces with constant conserved quantities and ask for a metric for measuring distances between isosurfaces.Because the Wasserstein distance measures distances between distributions rather than disjoint isosurfaces, it can easily generalize to noisy or partially observed data where the trajectories are no longer strictly disjoint.An alternative choice of metric, e.g. the Hausdorff distance between sets in a metric space [46], does not have this nice property and would be much more susceptible to noise or partial observations.On the other hand, because the Wasserstein distance distinguishes between different distributions, it is susceptible to other forms of corruption such as sampling inhomogeneity (discussed in Appendix G.3).In most cases, we believe that this tradeoff is worth it to retain the high degree of robustness as long as one is aware of the limitations.

Appendix F Comparison with Direct Fitting Methods
Only a few alternative approaches [14][15][16] have been proposed for identifying conservation laws from trajectory samples without time information (e.g.any methods that require time derivative estimates of the system state are not applicable).These alternatives all fall into a broad category that we term "direct fitting" methods, which attempt to directly fit a parameterized function of the system state to be constant over each trajectory.This can be accomplished using a variety of methods and optimization objectives but all generally require that the system state is fully observed and that the observed trajectories have low noise.These restrictions ensure that the basic assumption of direct fitting methods-that there exists a well-defined function from the observed state to the conserved quantities-is fulfilled.AI Poincaré [14] uses a symbolic regression approach to directly obtain an interpretable expression for the conserved quantities in terms of a library of symbolic expressions.While this method can sometimes give the exact expected conservation laws, AI Poincaré relies heavily on the applicability of symbolic regression, and, in cases where there is no simple symbolic expression, it will often fail.For example, for the planar gravitational dynamics dataset (named the "Kepler problem" in [14]), AI Poincaré fails to identify the conserved quantity associated to the orientation angle of the orbit due to its somewhat awkward symbolic representation [14].
The two remaining direct fitting methods, Siamese Neural Network (SNN) [15] and Con-servNet [16], both use a neural network as the parameterized function to fit the trajectory data.One important limitation of both of these approaches is their inability to discover more than a single conserved quantity per dataset.In both works [15,16], identifying a second conserved quantity requires generating a new dataset Table F2 Performance comparison of our method alongside deep learning-based direct fitting methods on the planar gravitational dynamics dataset.We show the rank correlation ρ of the embeddings with the ground truth conserved quantities.Note that the directing fitting approaches, ConservNet and Siamese Neural Network (SNN), are designed to discover a single conserved quantity and cannot handle multiple conserved quantities without dataset regeneration.Thus, with only a single dataset, we only expect the direct fitting methods to learn a single conserved quantity.For the direct fitting methods, when fitting to conserved quantities involving the orientation angle ϕ, we choose a constant ϕ 0 that gives the maximum ρ with conserved quantities of the form a cos(ϕ − ϕ 0 ).while holding the first discovered conserved quantity fixed.Generating such a dataset as part of the method pipeline assumes both access to and fine-grained control of the data generating process, which is often not available in practice.We benchmark our manifold learning approach against these two methods on the planar gravitational dynamics dataset, variations of which (under the name "Motion in a central potential" and "Kepler problem") were also previously studied in both works [15,16].Because we only use a single dataset, we only expect these direct fitting methods to identify a single conserved quantity.We also compare these methods on the noisy and partially observed versions of the dataset to understand their limitations.Implementations for both direct fitting methods were adapted from code released by [16] since reference code for [15] is unavailable.Following [16], our benchmark uses a simple fully connected neural network with four hidden layers of size 320, Mish activations [47], and an Adam optimizer [48] with learning rate 5 × 10 −5 .We use a batch size of 1000 for the SNN and 200 (fixed by the number of samples per trajectory) for ConservNet.We train each network for 1000 epochs4 and then compute the rank correlation ρ of the fitted function with the ground truth conserved quantities (Table F2).
The results show that, as expected, the direct fitting methods only learned a single conserved quantity for the dataset (Table F2).Even for the single conserved quantity identified by SNN and ConservNet, we see that our manifold learning method outperforms both direct fitting methods in all settings while also identifying all three conserved quantities.Training for SNN and Con-servNet took between 40-50 minutes for 1000 epochs on a single RTX 2080 Ti GPU and does not appear to benefit significantly from using larger batch sizes and additional compute resources.In comparison, our method-which is limited primarily by the Wasserstein distance estimation (Appendix G.1)-takes around 40 seconds to run on eight RTX 2080 Ti GPUs and 5-6 minutes on a single RTX 2080 Ti GPU.

Appendix G Additional
Discussion & Limitations

G.1 Time & Space Complexity
For a dataset with N trajectories and with S samples per trajectory, our approach is currently limited by the computational cost of estimating the 2-Wasserstein distance for all pairs of trajectories (Table G3), giving a time complexity of O(N 2 S 2 log S/ε 2 ) (Appendix B.1).To scale to much larger datasets, we have several choices to improve the run time performance.One simple adjustment to speed up the convergence of the Sinkhorn algorithm is to allow for a significantly larger target regularization parameter γ.The result, in fact, interpolates between the Wasserstein metric (γ = 0) and a maximum mean discrepancy (MMD) metric (γ = ∞) [35].However, to improve the time complexity of our approach as we scale to large S, we may have to consider more approximate approaches.A popular option is to first subsample the data to form minibatches of size s ≪ S, solve the much smaller optimal transport problem, and then average the resulting estimates for the Wasserstein distance [49,50].For a fixed number of epochs of averaging, this approach would give us linear scaling O(S) with the number of samples S.
To improve the scaling with the number of trajectories N , we may be able to take advantage of the sparse structure of the kernel matrix K when the total number of conserved quantities n ≪ N .That is, when the dimension of the embedded manifold (corresponding to the conserved quantities) is low, we expect each trajectory to have relatively few nearest neighbors, so most entries of the kernel matrix K will be very close to zero for an appropriately chosen Gaussian kernel.If that is the case, we can construct the kernel matrix K as a sparse matrix, e.g. by starting with a very coarse approximation for the pairwise Wasserstein distances and then obtaining finer estimates only for nearby trajectories or by constructing a k-nearest neighbor tree [51] (which takes O(kN log N ) time).With a sparse kernel matrix K, the diffusion map will only take O(N ) time and use O(N ) space.
In the future, by incorporating these approximations and algorithmic improvements, we expect to be able to adapt our approach to achieve linear scaling in time and space for very large datasets.

G.2 Sample Complexity
To understand the sample complexity of our approach, we need to consider the effect of both the number of trajectories N and the number of samples per trajectory S.
Table G3 Run times for computing the pairwise Wasserstein distances for each dataset.In addition to the run times, we also list the number of trajectories N in the dataset, the number of samples S per trajectory, and the dimension d of the phase space.We used eight RTX 2080 Ti GPUs for all computations.† Because we are interested in local conservation laws for the KdV equation dataset, the phase space is treated differently (see Sec.The number of samples S determines the accuracy of the estimated pairwise Wasserstein distances.For the Sinkhorn algorithm, the approximation error is [52] O S −1/2 e κ/γ 1 + γ −⌊d/2⌋ , (G11) where d is the dimension of the data, γ is the entropy regularization parameter, and κ is a datadependent constant.That is, the error in our Wasserstein distance estimates scales as 1/ √ S. Also, for γ ≥ 1, notice that the dimension d has no significant influence on the error bound.Furthermore, the relevant dimension d in the bound is likely to be some measure of the intrinsic dimension of the data [53], allowing us to obtain good Wasserstein distance estimates even for some systems with high dimensional phase spaces (Sec.4.5).
Assuming we have accurate Wasserstein distance estimates, the spectral error of the diffusion map is O(N −2/(8+n) ) for a manifold of dimension n (i.e. the number of conserved quantities) [54].Thus, for integrable PDE systems like the KdV equation (Sec.4.6) with an infinite number of conserved quantities, we must restrict ourselves to considering local conserved quantities to have a reasonable chance of reconstructing a useful embedding.

G.3 Robustness & Sampling Inhomogeneity
Because the Wasserstein distance is a metric over distributions, it has great robustness properties (Sec.E).However, the exact same qualities that provide this robustness also make the Wasserstein metric sensitive to sampling inhomogeneity between trajectories.That is, if two trajectories with the same conserved quantities are sampled in a way such that their distributions over the measured phase space differ, then the Wasserstein metric will treat them as different distributions, which could lead to spurious conserved quantities being identified by our method.Note that this does not include changes in the sampling process that are uniform across the trajectories and only a function of the phase space, e.g.sampling the state of the pendulum more often when it is near to the bottom of its swing, or sampling the position of a planet more often when it is farther from the sun for all trajectories.One example of the relevant kind of inhomogeneity is sampling trajectories for too short a time such that the physical measure (Eq.4) is not well approximated.Trajectories sampled over too short a time can essentially lead our method to believe that there are additional conservation laws preventing the system state from exploring areas of phase space that may in fact be reachable for a longer trajectory.On the hand, this can be considered a feature rather than a bug in the sense that, by choosing the time over which to sample our trajectories, we are providing our method with a time scale for our conserved quantities, allowing us to probe approximately conserved quantities that appear invariant over shorter time scales (see Appendix C for an example).
If we have a fully observed phase space and low noise, then this effect can be mitigated by choosing an alternative metric, such as a Hausdorff distance [46], that does not have a strong dependence on the sampling distribution.This would essentially take advantage of the fact that any overlap between two trajectories in a fully observed phase space implies that both have the same conserved quantities.However, for noisy data or a partially observed phase space, observed overlap between measured trajectories could also be due to noise or hidden variables.Again, the same qualities that make the Hausdorff distance robust to sampling inhomogeneity also make it very sensitive to noise and partial observations.
In other words, the Wasserstein metric is able to distinguish trajectories with different conserved quantities even in a noisy partially observed phase space precisely because it can sense differences in the distributions of the trajectories.

G.4 Physical Measures, Dissipation, & Unbounded Dynamics
The assumption that the dynamical system admits a physical measure (Sec.2.1.2) is a fairly natural one that allows us to characterize conserved quantities using invariant measures.In fact, for Hamiltonian dynamics that already admit a canonical Liouville measure [55], assuming that trajectories are ergodic on the isosurfaces of conserved quantities is essentially Boltzmann's famous ergodic hypothesis [56,57], which underpins much of statistical mechanics but is notoriously difficult to prove in general.
For dissipative systems, our assumption as stated is false since dissipation generally destroys ergodicity.However, dissipative systems still have attractors that admit physical measures [24,58].Therefore, for a dissipative system, rather than characterizing ergodic measures on isosurfaces, our method would instead be characterizing the various attractors of the dynamics, including fixed points, limit cycles, and chaotic attractors.Our experiment with the oscillating Turing pattern (Sec.4.5) is precisely such a dissipative system with a continuous set of chaotic attractors parameterized by a spatial phase angle η.
One clear example where our assumption of a physical measure is violated in a meaningful way is when the dynamics are unbounded.For example, hyperbolic orbits following planar gravitation dynamics escape to infinity and do not converge to any invariant measure.In this case, naively applying our approach would not lead to meaningful results since, for any fixed sampling time, two different trajectories from the same hyperbolic orbit will sample very different portions of phase space depending on initial conditions.One potential solution for this issue of unbounded dynamics is to first compactify the phase space.
For example, if we compactify the traditional four dimensional Euclidean phase space of planar gravitational dynamics into a four dimensional real projective space with an attached metric (e.g. the Fubini-Study metric) [59], we have effectively made the dynamics bounded and therefore amenable to our approach.Moving to this projective representation, the physical measures on the elliptical orbits would remain largely unchanged while trajectories on hyperbolic orbits would approach a fixed point at the "line at infinity" (which would be a finite "distance" away due to the choice of metric on this compactified space).Our method would then be able to characterize the elliptical orbits by the usual three conserved quantities and hyperbolic orbits by their limiting long time behavior corresponding to their final velocity vectors as they escape to infinity.

G.5 Trajectory Diversity & Correlated Conserved Quantities
One fundamental limitation, which affects not only our method but also other approaches for discovering conservation laws [12][13][14][15][16], is the inability to distinguish between constraints and correlated conserved quantities.In general, constraints on the system dynamics are constant across all possible initial conditions, whereas conserved quantities vary based on initial condition.However, if there is a lack of diversity in the initial conditions for the trajectories in our dataset, then it is possible to mistake two highly correlated conserved quantities for an additional constraint on the system.For example, if the planar gravitational dynamics dataset only contained nearly circular orbits, then the angular momentum and the energy of the orbits will become highly correlated and the orientation angle will be essentially meaningless.Given such a dataset of circular trajectories and no additional information, it is not possible to distinguish between a planar gravitational dynamics dataset with poor trajectory diversity and a dataset from a system that is constrained to circular orbits with a single conserved quantity (associated with the radius of the orbits).
As such, any general method for discovering conserved quantities will treat two highly correlated conserved quantities as a single conserved quantity.
Fig.1Proposed non-parametric method for discovering conservation laws illustrated using a simple pendulum example.(a) First, we collect and normalize the trajectory data from the dynamical system.Two example trajectories are highlighted in red and blue.(b) Then, we use the Wasserstein metric from optimal transport to compute the distance between each pair of trajectories and construct a distance matrix.For the two example trajectories, the optimal transport plan is shown as lines connecting pairs of points.The constructed distance matrix is plotted with color representing the computed Wasserstein distance between each pair of trajectories.The computed distance between the two example trajectories is marked (black dots) on the distance matrix plot.(c) An embedding of the shape space manifold C is extracted from the distance matrix using diffusion maps.The embedding plot is colored by the conserved energy of the pendulum E. The points corresponding to the two example trajectories are marked in red and blue.(d) Finally, a heuristic score (Appendix A) is used to select relevant components.In this case, only component 1 is relevant, corresponding to a single conserved quantity-the energy E. Again, the embedding plot is colored by E, and the two example trajectories are marked in red and blue.

Fig. 2
Fig. 2 Identifying the conserved energy for the simple harmonic oscillator (SHO).(a) The SHO has two degrees of freedom: position q and momentum p.(b) Sample trajectories from the SHO dataset show sample points plotted in the 2D phase space (q, p).(c) The heuristic score (with cutoff 0.6) correctly identifies that the first embedding component extracted by the diffusion map is the only relevant component.(d) The extracted first component closely matches the analytically predicted first component (Eq.15) for the SHO (R 2 = 0.9995).(e) Next, consider the SHO dataset with a partially observed phase space containing position only.(f) For each sample trajectory, the sample points are shown as a histogram.(g) The heuristic score is still able to identify the first component as relevant, and (h) this first component matches the analytic prediction (R 2 = 0.9961).

5 2 )Fig. 3
Fig. 3 Identifying the conserved energy for the simple pendulum.(a) The simple pendulum has two degrees of freedom: angular position θ and angular velocity ω.(b) Sample trajectories show sample points plotted in the 2D phase space (θ, ω).(c) The heuristic score (with cutoff 0.6) correctly identifies that the first embedding component extracted by the diffusion map is the only relevant component, and (d) the extracted first component is monotonically related to the energy (rank correlation ρ = 0.9997).(e, f) With the addition of σ = 0.5 Gaussian noise to simulate measurement noise, (g) the heuristic score is still able to identify the first component as relevant, and (h) this first component corresponds well to the energy (ρ = 0.9978).

Fig. 4
Fig. 4 Identifying conserved quantities for planar gravitational dynamics.(a) Planar gravitational dynamics has four degrees of freedom: position vector r and momentum vector p.(b) Sample trajectories show sample points plotted in 2D slices of the 4D phase space consisting of position r and momentum p. (c) The heuristic score (with cutoff 0.6) identifies three relevant embedding components corresponding to the three independent conserved quantities.(d, e) Components 1 and 2 embed the semimajor axis vector a with magnitude a = −1/2E related to the energy and orientation given by the angle ϕ.(f) Component 6 corresponds to the angular momentum L.

2 Fig. 5
Fig. 5 Identifying conserved quantities for the double pendulum.(a) The double pendulum has four degrees of freedom: angular positions θ 1 , θ 2 and angular velocities ω 1 , ω 2 .(b) Sample trajectories show sample points plotted in 2D slices of the 4D phase space.(c) The heuristic score (with cutoff 0.6) identifies one relevant embedding component corresponding to (e) the total energy E. (d) However, if we restrict the embedding to trajectories with first component v 1 < −1 (i.e.low energy trajectories) and renormalize the embedding, we find (f-h) two conserved quantities corresponding to the energies E ± of the two decoupled low energy modes.The gray points in Figures 5f-h correspond to the high energy trajectories (first component v 1 > −1) which are not relevant when considering the low energy non-chaotic phase of the double pendulum.

Fig. 6
Fig. 6 Identifying the conserved spatial phase for the oscillating Turing pattern system.(a) An example trajectory, with randomly sampled states u(x) and v(x) plotted, illustrates the high dimensional nature of the problem.(b)The heuristic score (with cutoff 0.6) identifies two relevant components, but on further examination, (c) we see that there is just a single conserved angle, corresponding to the spatial phase η of the Turing pattern, that needs to be embedded in two dimensions due to its topology.

Fig. 7
Fig. 7 Identifying three local conserved quantities of the Korteweg-De Vries (KdV) equation.(a) An example trajectory from the KdV dataset shows the high dimensional raw sampled states u(x).(b) To focus on local conserved quantities, we extract a distribution of the local features u(x), ∆u(x) from the raw states, removing the explicit spatial label.The plot shows the local feature distributions for a few sample states.(c) The heuristic score (with cutoff 0.6) correctly identifies three relevant components corresponding to (d-f) the three local conserved quantities (Eq.24).
Fig. A1 Breakdown of the heuristic score and illustration of redundant embedding components from the planar gravitational dynamics experiment.(a)The relative length scale l i /l 1 for each embedding component is computed from the corresponding eigenvalue λ i of the Laplacian matrix (Eq.A1).(b) The unpredictability measure m i for each component is computed using a nearest neighbor estimator[27].(c) The combined score m i l i /l 1 is the product of the relative length scale and the unpredictability measure.(d-f) The components 3, 4, and 5 are identified by the unpredictability measure as redundant.If we examine these three components, we find that they together embed a second order angular mode of components 1 and 2 (Figs.4d,e).In particular, the embedding is shaped like the surface of a cone with the height (or radial distance) roughly corresponding to the semi-major axis a and the angle around the cone corresponding to ϕ mod π, a second order mode of the orientation angle ϕ.

Fig.
Fig.A2Example from the planar gravitational dynamics experiment of identifying the number of conserved quantities (i.e. the embedding size).(a) Sweeping the cutoff value from 0 to 1, we find plateaus indicating robustness at embedding size 2 and 3. Note that there is a spurious plateau at the maximum embedding size 20.(b) A histogram of the embedding sizes confirms that the number of conserved quantities is likely to be 3.

Fig
Fig.B3Identifying conserved quantities for planar gravitational dynamics using Isomap.(a) When applying Isomap, the effective length scale λ Isomap can be used to estimate manifold dimensionality but it only clearly identifies two out of the three conserved quantities.(b, c) Similarly to the diffusion map embedding, Isomap components 1 and 2 embed the semi-major axis vector a with magnitude a = −1/2E related to the energy and orientation given by the angle ϕ.These components have high rank correlation ρ = 0.994 (ρ = 0.992) with a cos ϕ (a sin ϕ).(d) Component 3 roughly corresponds to the angular momentum L but is noisier than the embedding identified by diffusion maps (Fig.4f) and has a much lower rank correlation ρ = 0.723 with L.
Fig. C5 Effect of varying γ and τ on the identification of an approximate conservation law for the Langevin harmonic oscillator.Varying γ, τ , and the sampling time T , we show (a) the rank correlation ρ of the first embedding component with the time-averaged energy ⟨E⟩ and (b) the prominence s 1 − max({s i } i>1 ) of the score of the first component s 1 over the highest score of the remaining components s i .

2 Fig.
Fig. D6 Nonlinear periodic orbit of the double pendulum.(a) The four red highlighted points in the extracted embedding correspond to (b, c) a periodic orbit of the double pendulum that is connected to but well outside of the linear coupled oscillator regime.
Simple Pendulum: Position Only (b) Simple Pendulum: Position Only + Noise (c) Planar Gravitational Dynamics: Position Only (d) Planar Gravitational Dynamics: Noise

Fig.
Fig.E7Additional experiments illustrating the robustness of our approach.(a) For the simple pendulum system, even when provided only angle θ measurement data (without angular velocity ω), our method is able to identify a single relevant component corresponding to the energy the pendulum (ρ = 0.998).(b) If we then also add σ = 0.5 Gaussian noise, we can still achieve a similar result (ρ = 0.996).For planar gravitational dynamics, our method also performs well given (c) only position r data or (d) with σ = 0.5 Gaussian noise, correctly identifying the three conserved quantities.

Fig
Fig. E8 Identifying conserved quantities for the double pendulum from only position data.(a) The heuristic score (with cutoff 0.6) identifies one relevant embedding component corresponding to (b) the total energy E. (c) However, if we restrict the embedding to trajectories with first component v 1 < −1 (i.e.low energy trajectories) and renormalize the embedding, we find (d-f) two conserved quantities corresponding to the energies E ± of the two decoupled low energy modes.The gray points in Figures E8d-f correspond to the high energy trajectories (first component v 1 > −1) which are not relevant when considering the low energy non-chaotic phase of the double pendulum.
, Andrew Ma, Samuel Kim, Charlotte Loh, and Ruba Houssami.P.Y.L. gratefully acknowledges the support of the Eric and Wendy Schmidt AI in Science Postdoctoral Fellowship, a Schmidt Futures program.This research is supported in part by the U.S. Department of Defense through the National Defense Science & Engineering Graduate Fellowship Program; the National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, http://iaifi.org/);the U.S. Army Research Office through the Institute for Soldier Nanotechnologies at MIT under Collaborative Agreement Number W911NF-18-2-0048; the Air Force Office of Scientific Research under the award number FA9550-21-1-0317; and the United States Air Force Research Laboratory and the United States Air Force Artificial Intelligence Accelerator under Cooperative Agreement Number FA8750-19-2-1000.The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the United States Air Force or the U.S. Government.The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.
Table E1 Rank correlations ρ of linear fits with ground truth conserved quantities for the additional experiments.* The rank correlations for the low energy approximately conserved mode energies E ± are computed on the restricted set of trajectories with first embedding component v 1 < −1.
Bolded numbers indicate the best rank correlations for each conserved quantity. 4.6).