Efficiently Parameterized Neural Metriplectic Systems

## Abstract

Abstract Metriplectic systems are learned from data in a way that scales quadratically in both the size of the state and the rank of the metriplectic data. Besides being provably energy conserving and entropy stable, the proposed approach comes with approximation results demonstrating its ability to accurately learn metriplectic dynamics from data as well as an error estimate indicating its potential for generalization to unseen timescales when approximation error is low. Examples are provided which illustrate performance in the presence of both full state information as well as when entropic variables are unknown, confirming that the proposed approach exhibits superior accuracy and scalability without compromising on model expressivity.

## 1 Introduction

The theory of metriplectic, also called GENERIC, systems [, ] provides a principled formalism for encoding dissipative dynamics in terms of complete thermodynamical systems that conserve energy and produce entropy. By formally expressing the reversible and irreversible parts of state evolution with separate algebraic brackets, the metriplectic formalism provides a general mechanism for maintaining essential conservation laws while simultaneously respecting dissipative effects. Thermodynamic completeness implies that any dissipation is caught within a metriplectic system through the generation of entropy, allowing for a holistic treatment which has already found use in modeling fluids [, ], plasmas [, ], and kinetic theories [, ].

From a machine learning point of view, the formal separation of conservative and dissipative effects makes metriplectic systems highly appealing for the development of phenomenological models. Given data which is physics-constrained or exhibits some believed properties, a metriplectic system can be learned to exhibit the same properties with clearly identifiable conservative and dissipative parts. This allows for a more nuanced understanding of the governing dynamics via an evolution equation which reduces to an idealized Hamiltonian system as the dissipation is taken to zero. Moreover, elements in the kernel of the learned conservative part are immediately understood as Casimir invariants, which are special conservation laws inherent to the phase space of solutions, and are often useful for understanding and exerting control on low-dimensional structure in the system. On the other hand, the same benefit of metriplectic structure as a “direct sum” of reversible and irreversible parts makes it challenging to parameterize in an efficient way, since delicate degeneracy conditions must be enforced in the system for all time. In fact, there are no methods at present for learning general metriplectic systems which scale optimally with both dimension and the rank of metriplectic data—an issue which this work directly addresses.

Precisely, metriplectic dynamics on the finite or infinite dimensional phase space $\mathcal{P}$ are generated by a free energy function(al) $F:\mathcal{P}\to\mathbb{R}$, $F=E+S$ defined in terms of a pair $E,S:\mathcal{P}\to\mathbb{R}$ representing energy and entropy, respectively, along with two algebraic brackets $\mathopen{}\mathclose{{}\left\{\cdot,\cdot}\right\},\mathopen{}\mathclose{{}
\left[\cdot,\cdot}\right]:C^{\infty}\mathopen{}\mathclose{{}\left(\mathcal{P}}
\right)\times C^{\infty}\mathopen{}\mathclose{{}\left(\mathcal{P}}\right)\to C
^{\infty}\mathopen{}\mathclose{{}\left(\mathcal{P}}\right)$ which are bilinear derivations on $C^{\infty}\mathopen{}\mathclose{{}\left(\mathcal{P}}\right)$ with prescribed symmetries and degeneracies $\mathopen{}\mathclose{{}\left\{S,\cdot}\right\}=\mathopen{}\mathclose{{}\left[
E,\cdot}\right]=0$. Here $\mathopen{}\mathclose{{}\left\{\cdot,\cdot}\right\}$ is an antisymmetric Poisson bracket, which is a Lie algebra realization on functions, and $\mathopen{}\mathclose{{}\left[\cdot,\cdot}\right]$ is a degenerate metric bracket which is symmetric and positive semi-definite. When $\mathcal{P}\subset\mathbb{R}^{n}$ for some $n>0$, these brackets can be identified with symmetric matrix fields $~{}L:\mathcal{P}\to\mathrm{Skew}_{n}\mathopen{}\mathclose{{}\left(\mathbb{R}}
\right),~{}M:\mathcal{P}\to\mathrm{Sym}_{n}\mathopen{}\mathclose{{}\left(
\mathbb{R}}\right)$ satisfying $\mathopen{}\mathclose{{}\left\{F,G}\right\}=\nabla F\cdot~{}L\nabla G$ and $\mathopen{}\mathclose{{}\left[F,G}\right]=\nabla F\cdot~{}M\nabla G$
for all functions $F,G\in C^{\infty}\mathopen{}\mathclose{{}\left(\mathcal{P}}\right)$ and all states $~{}x\in\mathcal{P}$. Using the degeneracy conditions along with $\nabla~{}x=~{}I$ and abusing notation slightly
then leads the standard equations for metriplectic dynamics,

$$ $\displaystyle\dot{~{}x}$ $\displaystyle=\mathopen{}\mathclose{{}\left\{~{}x,F}\right\}+\mathopen{} \mathclose{{}\left[~{}x,F}\right]=\mathopen{}\mathclose{{}\left\{~{}x,E}\right \}+\mathopen{}\mathclose{{}\left[~{}x,S}\right]=~{}L\nabla E+~{}M\nabla S,$ $$

which are provably energy conserving and entropy producing. To see why this is the case, recall that $~{}L^{\intercal}=-~{}L$. It follows that the infinitesimal change in energy satisfies

$$ $\displaystyle\dot{E}$ $\displaystyle=\dot{~{}x}\cdot\nabla E=~{}L\nabla E\cdot\nabla E+~{}M\nabla S \cdot\nabla E=-~{}L\nabla E\cdot\nabla E+\nabla S\cdot~{}M\nabla E=0,$ $$

and therefore energy is conserved along the trajectory of $~{}x$. Similarly, the fact that $~{}M^{\intercal}=~{}M$ is positive semi-definite implies that

$$ $\displaystyle\dot{S}$ $\displaystyle=\dot{~{}x}\cdot\nabla S=L\nabla E\cdot\nabla S+~{}M\nabla S\cdot \nabla S=-\nabla E\cdot~{}L\nabla S+~{}M\nabla S\cdot\nabla S=\mathopen{} \mathclose{{}\left|\nabla S}\right|_{~{}M}^{2}\geq 0,$ $$

so that entropy is nondecreasing along $~{}x$ as well. Geometrically, this means that the motion of a trajectory $~{}x$ is everywhere tangent to the level sets of energy and transverse to those of entropy, reflecting the fact that metriplectic dynamics are a combination of noncanonical Hamiltonian ($~{}M=~{}0$) and generalized gradient ($~{}L=~{}0$) motions. Note that these considerations also imply the Lyapunov stability of metriplectic trajectories, as can be seen by taking $E$ as a Lyapunov function. Importantly, this also implies that metriplectic trajectories which start in the (often compact) set $K=\{~{}x\,|\,E(~{}x)\leq E(~{}x_{0})\}$ remain there for all time.

In phenomenological modeling, the entropy $S$ is typically chosen from Casimirs of the Poisson bracket generated by $~{}L$, i.e. those quantities $C\in C^{\infty}\mathopen{}\mathclose{{}\left(\mathcal{P}}\right)$ such that $~{}L\nabla C=~{}0$. On the other hand, the method which will be presented here, termed neural metriplectic systems (NMS), allows for all of the metriplectic data $~{}L,~{}M,E,S$ to be approximated simultaneously, removing the need for Casimir invariants to be known or assumed ahead of time. The only restriction inherent to NMS is that the metriplectic system being approximated is nondegenerate (c.f. Definition 3.1), a mild condition meaning that the gradients of energy and entropy cannot vanish at any point $~{}x\in\mathcal{P}$ in the phase space. It will be shown that NMS alleviates known issues with previous methods for metriplectic learning, leading to easier training, superior parametric efficiency, and better generalization performance.

#### Contributions.

The proposed NMS method for learning metriplectic models offers the following advantages over previous state-of-the-art: (1) It approximates arbitrary nondegenerate metriplectic dynamics with optimal quadratic scaling in both the problem dimension $n$ and the rank $r$ of the irreversible dynamics. (2) It produces realistic, thermodynamically consistent entropic dynamics from unobserved entropy data. (3) It admits universal approximation and error accumulation results given in Proposition 3.7 and Theorem 3.9. (4) It yields exact energy conservation and entropy stability by construction, allowing for effective generalization to unseen timescales.

## 2 Previous and Related Work

Previous attempts to learn metriplectic systems from data separate into “hard” and “soft” constrained methods. Hard constrained methods enforce metriplectic structure by construction, so that the defining properties of metriplecticity cannot be violated. Conversely, methods with soft constraints relax some aspects of metriplectic structure in order to produce a wider model class which is easier to parameterize. While hard constraints are the only way to truly guarantee appropriate generalization in the learned surrogate, the hope of soft constrained methods is that the resulting model is “close enough” to metriplectic that it will exhibit some of the favorable characteristics of metriplectic systems, such as energy and entropy stability. Some properties of the methods compared in this work are summarized in Table 3.3.

#### Soft constrained methods.

Attempts to learn metriplectic systems using soft constraints rely on relaxing the degeneracy conditions $~{}L\nabla S=~{}M\nabla E=~{}0$. This is the approach taken in [], termed SPNN, which learns an almost-metriplectic model parameterized with generic neural networks through a simple $L^{2}$ penalty term in the training loss, $\mathcal{L}_{\mathrm{pen}}=\mathopen{}\mathclose{{}\left|~{}L\nabla E}\right|^
{2}+\mathopen{}\mathclose{{}\left|~{}M\nabla S}\right|^{2}.$
This widens the space of allowable network parameterizations for the approximate operators $~{}L,~{}M$. While the resulting model violates the first and second laws of thermodynamics, the authors show that reasonable trajectories are still obtained, at least when applied within the range of timescales used for training. A similar approach is taken in [], which targets larger problems and develops an almost-metriplectic model reduction strategy based on the same core idea.

#### Hard constrained methods.

Perhaps the first example of learning metriplectic systems from data was given in [] in the context of system identification. Here, training data is assumed to come from a finite element simulation, so that the discrete gradients of energy and entropy can be approximated as $~{}\nabla E(~{}x)=~{}A~{}x,\nabla S(~{}x)=~{}B~{}x$. Assuming a fixed form for $~{}L$ produces a constrained learning problem for the constant matrices $~{}M,~{}A,~{}B$ which is solved to yield a provably metriplectic surrogate model. Similarly, the work [] learns $~{}M,E$ given $~{}L,S$ by considering a fixed block-wise decoupled form which trivializes the degeneracy conditions, i.e. $~{}L=[\star\,\,~{}0;~{}0\,\,~{}0]$ and $~{}M=[~{}0\,\,~{}0;~{}0\,\,\star]$.
This line of thought is continued in [] and [], both of which learn metriplectic systems with neural network parameterizations by assuming this decoupled block structure. A somewhat broader class of metriplectic systems are considered in [] using tools from exterior calculus, with the goal of learning metriplectic dynamics on graph data. This leads to a structure-preserving network surrogate which scales linearly in the size of the graph domain, but also cannot express arbitrary metriplectic dynamics due to the specific choices of model form for $~{}L,~{}M$.

A particularly inspirational method for learning general metriplectic systems was given in [] and termed GNODE, building on parameterizations of metriplectic operators developed in []. GNODE parameterizes learnable reversible and irreversible bracket generating matrices $~{}L,~{}M$ in terms of state-independent tensors $~{}\xi\in\mathopen{}\mathclose{{}\left(\mathbb{R}^{n}}\right)^{\otimes 3}$ and $~{}\zeta\in\mathopen{}\mathclose{{}\left(\mathbb{R}^{n}}\right)^{\otimes 4}$: for $1\leq\alpha,\beta,\gamma,\mu,\nu\leq n$, the authors choose $L_{\alpha\beta}(~{}x)=\sum_{\gamma}\xi_{\alpha\beta\gamma}\partial^{\gamma}S$ and $M_{\alpha\beta}(~{}x)=\sum_{\mu,\nu}\zeta_{\alpha\mu,\beta\nu}\partial^{\mu}E
\partial^{\nu}E,$
where $\partial^{\alpha}F=\partial F/\partial x_{\alpha}$, $~{}\xi$ is totally antisymmetric, and $~{}\zeta$ is symmetric between the pairs $(\alpha,\mu)$ and $(\beta,\nu)$ but antisymmetric within each of these pairs.
The key idea here is to exchange the problem of enforcing degeneracy conditions $~{}L\nabla E=~{}M\nabla S=~{}0$ in matrix fields $~{}L,~{}M$ with the problem of enforcing symmetry conditions in tensor fields $~{}\xi,~{}\zeta$, which is comparatively easier but comes at the expense of underdetermining the problem. In GNODE, enforcement of these symmetries is handled by a generic learnable 3-tensor $\tilde{~{}\xi}\in\mathopen{}\mathclose{{}\left(\mathbb{R}^{n}}\right)^{\otimes
3}$ along with learnable matrices $~{}D\in\mathrm{Sym}_{r}\mathopen{}\mathclose{{}\left(\mathbb{R}}\right)$ and $~{}\Lambda^{s}\in\mathrm{Skew}_{n}\mathopen{}\mathclose{{}\left(\mathbb{R}}\right)$ for $1\leq s\leq r\leq n$, leading to the final parameterizations $\xi_{\alpha\beta\gamma}=\frac{1}{3!}\mathopen{}\mathclose{{}\left(\tilde{\xi}_
{\alpha\beta\gamma}-\tilde{\xi}_{\alpha\gamma\beta}+\tilde{\xi}_{\beta\gamma
\alpha}-\tilde{\xi}_{\beta\alpha\gamma}+\tilde{\xi}_{\gamma\alpha\beta}-\tilde
{\xi}_{\gamma\beta\alpha}}\right)$ and $\zeta_{\alpha\mu,\beta\nu}=\sum_{s,t}\Lambda_{\alpha\mu}^{s}D_{st}\Lambda_{
\beta\nu}^{t}.$
Along with learnable energy and entropy functions $E,S$ parameterized by multi-layer perceptrons (MLPs), the data $~{}L,~{}M$ learned by GNODE
guarantees metriplectic structure in the surrogate model and leads to successful learning of metriplectic systems in some simple cases of interest. However, note that this is a highly redundant parameterization requiring ${n\choose 3}+r{n\choose 2}+{{r+1}\choose 2}+2$ learnable scalar functions, which exhibits $\mathcal{O}(n^{3}+rn^{2})$ scaling in the problem size because of the necessity to compute and store ${n\choose 3}$ entries of $~{}\xi$ and $r{n\choose 2}$ entries of $~{}\Lambda$. Additionally, the assumption of state-independence in the bracket generating tensors $~{}\xi,~{}\zeta$ is somewhat restrictive, limiting the class of problems to which GNODE can be applied.

A related approach to learning metriplectic dynamics with hard constraints was seen in [], which proposed a series of GFINN architectures depending on how much of the information $~{}L,~{}M,E,S$ is assumed to be known. In the case that $~{}L,~{}M$ are known, the GFINN energy and entropy are parameterized with scalar-valued functions $f\circ~{}P_{\mathrm{ker}~{}A}$ where $f:\mathbb{R}^{n}\to\mathbb{R}$ ($E$ or $S$) is learnable and $~{}P_{\mathrm{ker}~{}A}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is orthogonal projection onto the kernel of the (known) operator $~{}A$ ($~{}L$ or $~{}M$). It follows that the gradient $\nabla\mathopen{}\mathclose{{}\left(f\circ~{}P_{\mathrm{ker}~{}A}}\right)=~{}P
_{\mathrm{ker}~{}A}\nabla f\mathopen{}\mathclose{{}\left(P_{\mathrm{ker}~{}A}}\right)$ lies in the kernel of $~{}A$, so that the degeneracy conditions are guaranteed at the expense of constraining the model class of potential energies/entropies. Alternatively, in the case that all of $~{}L,~{}M,E,S$ are unknown, GFINNs use learnable scalar functions $f$ for $E,S$ parameterized by MLPs as well as two matrix fields $~{}Q^{E},~{}Q^{S}\in\mathbb{R}^{r\times n}$ with rows given by $~{}q^{f}_{s}=\mathopen{}\mathclose{{}\left(~{}S^{f}_{s}\nabla f}\right)^{\intercal}$ for learnable skew-symmetric matrices $~{}S^{f}_{s}$, $1\leq s\leq r$, $f=E,S$. Along with two triangular $(r\times r)$ matrix fields $~{}T_{~{}L},~{}T_{~{}M}$, this yields the parameterizations
$~{}L(~{}x)=~{}Q^{S}(~{}x)^{\intercal}\mathopen{}\mathclose{{}\left(~{}T_{~{}L}
(~{}x)^{\intercal}-~{}T_{~{}L}(~{}x)}\right)~{}Q^{S}(~{}x)$ and $~{}M(~{}x)=~{}Q^{E}(~{}x)^{\intercal}\mathopen{}\mathclose{{}\left(~{}T_{~{}M}
(~{}x)^{\intercal}~{}T_{~{}M}(~{}x)}\right)~{}Q^{E}(~{}x)$,
which necessarily satisfy the symmetries and degeneracy conditions required for metriplectic structure. GFINNs are shown to both increase expressivity over the GNODE method as well as decrease redundancy, since the need for an explicit order-3 tensor field is removed and the reversible and irreversible brackets are allowed to depend explicitly on the state $~{}x$. However, GFINNs still exhibit cubic scaling through the requirement of $rn(n-1)+r^{2}+2=\mathcal{O}\mathopen{}\mathclose{{}\left(rn^{2}}\right)$ learnable functions, which is well above the theoretical minimum required to express a general metriplectic system and leads to difficulties in training the resulting models.

#### Model reduction.

Finally, it is worth mentioning the closely related line of work involving model reduction for metriplectic systems, which began with []. As remarked there, preserving metriplecticity in reduced-order models (ROMs) exhibits many challenges due to its delicate requirements on the kernels of the involved operators. There are also hard and soft constrained approaches: the already mentioned [] aims to learn a close-to-metriplectic data-driven ROM by enforcing degeneracies by penalty, while [] directly enforces metriplectic structure in projection-based ROMs using exterior algebraic factorizations. The parameterizations of metriplectic data presented here are related to those presented in [], although NMS does not require access to nonzero components of $\nabla E,\nabla S$.

## 3 Formulation and Analysis

The proposed formulation of the metriplectic bracket-generating operators $~{}L,~{}M$ used by NMS is based on the idea of exploiting structure in the tensor fields $~{}\xi,~{}\zeta$ to reduce the necessary number of degrees of freedom. In particular, it will be shown that the degeneracy conditions $~{}L\nabla S=~{}M\nabla E=~{}0$ imply more than just symmetry constraints on these fields, and that taking these additional constraints into account allows for a more compact representation of metriplectic data. Following this, results are presented which show that the proposed formulation is universally approximating on nondegenerate systems (c.f. Definition 3.1) and admits a generalization error bound in time.

### 3.1 Exterior algebra

Developing these metriplectic expressions will require some basic facts from exterior algebra, of which more details can be found in, e.g., [, Chapter 19]. The basic objects in the exterior algebra $\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge}}{{
\scriptstyle\wedge}}V$ over the vector space $V$ are multivectors, which are formal linear combinations of totally antisymmetric tensors on $V$. More precisely, if $I(V)$ denotes the two-sided ideal of the free tensor algebra $T(V)$ generated by elements of the form $~{}v\otimes~{}v$ ($~{}v\in V$), then the exterior algebra is the quotient space $\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge}}{{
\scriptstyle\wedge}}V\simeq T(V)/I(V)$ equipped with the antisymmetric wedge product operation $\wedge$. This graded algebra is equipped with natural projection operators $P^{k}:\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge}}{{
\scriptstyle\wedge}}V\to\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{
\textstyle\wedge}}{{\scriptstyle\wedge}}^{k}V$ which map between the full exterior algebra and the $k^{\mathrm{th}}$ exterior power of $V$, denoted $\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge}}{{
\scriptstyle\wedge}}^{k}V$, whose elements are homogeneous $k$-vectors. More generally, given an $n$-manifold $M$ with tangent bundle $TM$, the exterior algebra $\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge}}{{
\scriptstyle\wedge}}\mathopen{}\mathclose{{}\left(TM}\right)$ is the algebra of multivector fields whose fiber over $x\in M$ is given by $\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge}}{{
\scriptstyle\wedge}}T_{x}M$.

For the present purposes, it will be useful to develop a correspondence between bivectors $\mathsf{B}\in\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge
}}{{\scriptstyle\wedge}}^{2}\mathopen{}\mathclose{{}\left(\mathbb{R}^{n}}\right)$ and skew-symmetric matrices $~{}B\in\mathrm{Skew}_{n}\mathopen{}\mathclose{{}\left(\mathbb{R}}\right)$, which follows directly from Riesz representation in terms of the usual Euclidean dot product. More precisely, supposing that $~{}e_{1},...,~{}e_{n}$ are the standard basis vectors for $\mathbb{R}^{n}$, any bivector $\mathsf{B}\in\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge
}}{{\scriptstyle\wedge}}^{2}T\mathbb{R}^{n}$ can be represented as $\mathsf{B}=\sum_{i<j}B^{ij}~{}e_{i}\wedge~{}e_{j}$ where $B^{ij}\in\mathbb{R}$ denote the components of $\mathsf{B}$. Define a grade-lowering action of bivectors on vectors through right contraction (see e.g. Section 3.4 of []), expressed for any vector $~{}v$ and basis bivector $~{}e_{i}\wedge~{}e_{j}$ as
$\mathopen{}\mathclose{{}\left(~{}e_{i}\wedge~{}e_{j}}\right)\cdot~{}v=
\mathopen{}\mathclose{{}\left(~{}e_{j}\cdot~{}v}\right)~{}e_{i}-\mathopen{}
\mathclose{{}\left(~{}e_{i}\cdot~{}v}\right)~{}e_{j}.$ It follows that the action of $\mathsf{B}$ is equivalent to

$$ $\displaystyle\mathsf{B}\cdot~{}v$ $\displaystyle=\sum_{i<j}B^{ij}\mathopen{}\mathclose{{}\left(\mathopen{} \mathclose{{}\left(~{}e_{j}\cdot~{}v}\right)~{}e_{i}-\mathopen{}\mathclose{{} \left(~{}e_{i}\cdot~{}v}\right)~{}e_{j}}\right)=\sum_{i<j}B^{ij}v_{j}~{}e_{i}- \sum_{j<i}B^{ji}v_{j}~{}e_{i}=\sum_{i,j}B^{ij}v_{j}~{}e_{i}=~{}B~{}v,$ $$

where $~{}B^{\intercal}=-~{}B\in\mathbb{R}^{n\times n}$ is a skew-symmetric matrix representing $\mathsf{B}$, and we have re-indexed under the second sum and applied that $B^{ij}=-B^{ji}$ for all $i,j$. Since the kernel of this action is the zero bivector, it is straightforward to check that this string of equalities defines an isomorphism $\mathcal{M}:\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle\wedge}
}{{\scriptstyle\wedge}}^{2}\mathbb{R}^{n}\to\mathrm{Skew}_{n}\mathopen{}
\mathclose{{}\left(\mathbb{R}}\right)$ from the $2^{\mathrm{nd}}$ exterior power of $\mathbb{R}^{n}$ to the space of skew-symmetric $(n\times n)$-matrices over $\mathbb{R}$: in what follows, we will write $~{}B\simeq\mathsf{B}$ rather than $~{}B=\mathcal{M}(\mathsf{B})$ for notational convenience. Note that a correspondence in the more general case of bivector/matrix fields follows in the usual way via the fiber-wise extension of $\mathcal{M}$.

### 3.2 Learnable metriplectic operators

It is now possible to explain the proposed NMS formulation. First, note the following key definition which prevents the consideration of unphysical examples.

###### Definition 3.1 .

A metriplectic system on $K\subset\mathbb{R}^{n}$ generated by the data $~{}L,~{}M,E,S$ will be called nondegenerate provided $\nabla E,\nabla S\neq~{}0$ for all $~{}x\in K$.

With this, the NMS parameterizations for metriplectic operators are predicated on an algebraic result proven in Appendix A.

###### Lemma 3.2 .

Let $K\subset\mathbb{R}^{n}$. For all $~{}x\in K$, the operator $~{}L:K\to\mathbb{R}^{n\times n}$ satisfies $~{}L^{\intercal}=-~{}L$ and $~{}L\nabla S=~{}0$ for some $S:K\to\mathbb{R}$, $\nabla S\neq~{}0$, provided there exists a non-unique bivector field $\mathsf{A}:U\to\mathchoice{{\textstyle\bigwedge}}{{\bigwedge}}{{\textstyle
\wedge}}{{\scriptstyle\wedge}}^{2}\mathbb{R}^{n}$ and equivalent matrix field $~{}A\simeq\mathsf{A}$ such that

$$ $\displaystyle~{}L$ $\displaystyle\simeq\mathopen{}\mathclose{{}\left(\mathsf{A}\wedge\frac{\nabla S }{\mathopen{}\mathclose{{}\left|\nabla S}\right|^{2}}}\right)\cdot\nabla S= \mathsf{A}-\frac{1}{\mathopen{}\mathclose{{}\left|\nabla S}\right|^{2}}~{}A \nabla S\wedge\nabla S.$ $$

Similarly, for all $~{}x\in K$ a positive semi-definite operator $~{}M:K\to\mathbb{R}^{n\times n}$ satisfies $~{}M^{\intercal}=~{}M$ and $~{}M\nabla E=~{}0$ for some $E:K\to\mathbb{R}$, $\nabla E\neq~{}0$, provided there exists a non-unique matrix-valued $~{}B:K\to\mathbb{R}^{n\times r}$ and symmetric matrix-valued $~{}D:K\to\mathbb{R}^{r\times r}$ such that $r\leq n$ and

$$ $\displaystyle~{}M$ $\displaystyle=\sum_{s,t}D_{st}\mathopen{}\mathclose{{}\left(~{}b^{s}\wedge \frac{\nabla E}{\mathopen{}\mathclose{{}\left|\nabla E}\right|^{2}}}\right) \cdot\nabla E\,\otimes\,\mathopen{}\mathclose{{}\left(~{}b^{t}\wedge\frac{ \nabla E}{\mathopen{}\mathclose{{}\left|\nabla E}\right|^{2}}}\right)\cdot\nabla E$ $\displaystyle=\sum_{s,t}D_{st}\mathopen{}\mathclose{{}\left(~{}b^{s}-\frac{~{} b^{s}\cdot\nabla E}{\mathopen{}\mathclose{{}\left|\nabla E}\right|^{2}}\nabla E }\right)\mathopen{}\mathclose{{}\left(~{}b^{t}-\frac{~{}b^{t}\cdot\nabla E}{ \mathopen{}\mathclose{{}\left|\nabla E}\right|^{2}}\nabla E}\right)^{\intercal},$ $$

where $~{}b^{s}$ denotes the $s^{\mathrm{th}}$ column of $~{}B$. Moreover, using $~{}P^{\perp}_{f}=\mathopen{}\mathclose{{}\left(~{}I-\frac{\nabla f\,\nabla f^{
\intercal}}{\mathopen{}\mathclose{{}\left|\nabla f}\right|^{2}}}\right)$ to denote the orthogonal projector onto $\mathrm{Span}\mathopen{}\mathclose{{}\left(\nabla f}\right)^{\perp}$, these parameterizations of $~{}L,~{}M$ are equivalent to the matricized expressions $~{}L=~{}P^{\perp}_{S}~{}A~{}P^{\perp}_{S}$ and $~{}M=~{}P^{\perp}_{E}~{}B~{}D~{}B^{\intercal}~{}P^{\perp}_{E}.$

###### Remark 3.3 .

Observe that the projections appearing in these expressions are the minimum necessary for guaranteeing the symmetries and degeneracy conditions necessary for metriplectic structure. In particular, conjugation by $~{}P^{\perp}_{f}$ respects symmetry and ensures that both the input and output to the conjugated matrix field lie in $\mathrm{Span}(\nabla f)^{\perp}$.

Lemma 3.2 demonstrates specific parameterizations for $~{}L,~{}M$ that hold for any nondegenerate metriplectic data and are core to the NMS method for learning metriplectic dynamics. While generally underdetermined, these expressions are in a sense maximally specific given no additional information, since there is nothing available in the general metriplectic formalism to determine the matrix fields $~{}A,~{}B~{}D~{}B^{\intercal}$ on $\mathrm{Span}\mathopen{}\mathclose{{}\left(\nabla S}\right),\mathrm{Span}
\mathopen{}\mathclose{{}\left(\nabla E}\right)$, respectively. The following Theorem, also proven in Appendix A, provides a rigorous correspondence between metriplectic systems and these particular parameterizations.

###### Theorem 3.4 .

The data $~{}L,~{}M,E,S$ form a nondegenerate metriplectic system in the state variable $~{}x\in K$ if and only if there exist a skew-symmetric $~{}A:K\to\mathrm{Skew}_{n}\mathopen{}\mathclose{{}\left(\mathbb{R}}\right)$, symmetric postive semidefinite $~{}D:K\to\mathrm{Sym}_{r}\mathopen{}\mathclose{{}\left(\mathbb{R}}\right)$, and generic $~{}B:K\to\mathbb{R}^{n\times r}$ such that

$$ $\displaystyle\dot{~{}x}$ $\displaystyle=~{}L\nabla E+~{}M\nabla S=~{}P^{\perp}_{S}~{}A~{}P^{\perp}_{S} \nabla E+~{}P^{\perp}_{E}~{}B~{}D~{}B^{\intercal}~{}P^{\perp}_{E}\nabla S.$ $$

###### Remark 3.5 .

Note that the proposed parameterizations for $~{}L,~{}M$ are not one-to-one but properly contain the set of valid nondegenerate metriplectic systems in $E,S$, since the Jacobi identity on $~{}L$ necessary for a true Poisson manifold structure is not strictly enforced. For $1\leq i,j,k\leq n$, the Jacobi identity is given in components as $\sum_{\ell}L_{i\ell}\partial^{\ell}L_{jk}+L_{j\ell}\partial^{\ell}L_{ki}+L_{k
\ell}\partial^{\ell}L_{ij}=0.$
However, this requirement is not often enforced in algorithms for learning general metriplectic (or even symplectic) systems, since it is considered subordinate to energy conservation and it is well known that both qualities cannot hold simultaneously in general [].

###### Definition 3.1 .

###### Lemma 3.2 .

###### Remark 3.3 .

###### Theorem 3.4 .

###### Remark 3.5 .

### 3.3 Specific parameterizations

Now that Theorem 3.4 has provided a model class which is rich enough to express any desired metriplectic system, it remains to discuss what NMS actually learns. First, note that it is unlikely to be the case that $E,S$ are known a priori, so it is beneficial to allow these functions to be learnable alongside the governing operators $~{}L,~{}M$. For simplicity, energy and entropy $E,S$ are parameterized as scalar-valued MLPs with $\tanh$ activation, although any desired architecture could be chosen for this task. The skew-symmetric matrix field $~{}A=~{}A_{\mathrm{tri}}-~{}A_{\mathrm{tri}}^{\intercal}$ used to build $~{}L$ is parameterized through its strictly lower-triangular part $~{}A_{\mathrm{tri}}$ using a vector-valued MLP with output dimension ${n\choose 2}$, which guarantees that the mapping $~{}A_{\mathrm{tri}}\mapsto~{}A$ above is bijective.
Similarly, the symmetric matrix field $~{}D=~{}K_{\mathrm{chol}}~{}K_{\mathrm{chol}}^{\intercal}$ is parameterized through its lower-triangular Cholesky factor $~{}K_{\mathrm{chol}}$, which is a vector-valued MLP with output dimension ${r+1\choose 2}$. While this choice does not yield a bijective mapping $~{}K_{\mathrm{chol}}\mapsto~{}D$ unless, e.g., $~{}D$ is assumed to be positive definite with diagonal entries of fixed sign, this does not hinder the method in practice. In fact, $~{}D$ should not be positive definite in general, as this is equivalent to claiming that $~{}M$ is positive definite on vectors tangent to the level sets of $E$. Finally, the generic matrix field $~{}B$ is parameterized as a vector-valued MLP with output dimension $nr$. Remarkably, the exterior algebraic expressions in Lemma 3.2 require less redundant operations than the corresponding matricized expressions from Theorem 3.4, and therefore the expressions from Lemma 3.2 are used when implementing NMS. Figure 3.3 summarizes this information.

###### Remark 3.6 .

With these choices, the NMS parameterization of metriplectic systems requires only $(1/2)\mathopen{}\mathclose{{}\left((n+r)^{2}-(n-r)}\right)+2$ learnable scalar functions, in contrast to ${n\choose 3}+r{n\choose 2}+{r+1\choose 2}+2$ for the GNODE approach in [] and $rn(n-1)+r^{2}+2$ for the GFINN approach in []. In particular, NMS is quadratic in both $n,r$ with no decrease in model expressivity, in contrast to the cubic scaling of previous methods.

tableProperties of the metriplectic architectures compared.
Name
Physics Bias
Restrictive
Scale
NODE
None
No
Linear
SPNN
Soft
No
Quadratic
GNODE
Hard
Yes
Cubic
GFINN
Hard
No
Cubic
NMS
Hard
No
Quadratic

figureA visual depiction of the NMS architecture.

###### Remark 3.6 .

### 3.4 Approximation and error

Besides offering a compact parameterization of metriplectic dynamics, the expressions used in NMS also exhibit desirable approximation properties which guarantee a reasonable bound on state error over time. To state this precisely, first note the following universal approximation result proven in Appendix A.

###### Proposition 3.7 .

Let $K\subset\mathbb{R}^{n}$ be compact and $E,S:K\to\mathbb{R}$ be continuous such that $~{}L\nabla S=~{}M\nabla E=~{}0$ and $\nabla E,\nabla S\neq~{}0$ for all $~{}x\in K$. For any $\varepsilon>0$, there exist two-layer neural network functions $\tilde{E},\tilde{S}:K\to\mathbb{R},\tilde{~{}L}:K\to\mathrm{Skew}_{n}\mathopen
{}\mathclose{{}\left(\mathbb{R}}\right)$ and $\tilde{~{}M}:K\to\mathrm{Sym}_{n}\mathopen{}\mathclose{{}\left(\mathbb{R}}\right)$ such that $\nabla\tilde{E},\nabla\tilde{S}\neq~{}0$ on $K$, $\tilde{~{}M}$ is positive semi-definite, $\tilde{~{}L}\nabla\tilde{S}=\tilde{~{}M}\nabla\tilde{E}=~{}0$ for all $~{}x\in K$, and each approximate function is $\varepsilon$-close to its given counterpart on $K$. Moreover, if $~{}L,~{}M$ have $k\geq 0$ continuous derivatives on $K$ then so do $\tilde{~{}L},\tilde{~{}M}$.

###### Remark 3.8 .

The assumption $~{}x\in K$ of the state remaining in a compact set $V$ is not restrictive when at least one of $E,-S:\mathbb{R}^{n}\to\mathbb{R}$, say $E$, has bounded sublevel sets. In this case, letting $~{}x_{0}=~{}x(0)$ it follows from $\dot{E}\leq 0$ that
$E\mathopen{}\mathclose{{}\left(~{}x(t)}\right)=E\mathopen{}\mathclose{{}\left(
~{}x_{0}}\right)+\int_{0}^{t}\dot{E}\mathopen{}\mathclose{{}\left(~{}x(\tau)}
\right)\,d\tau\leq E\mathopen{}\mathclose{{}\left(~{}x_{0}}\right),$
so that the entire trajectory $~{}x(t)$ lies in the (closed and bounded) compact set $K=\{~{}x\,|\,E(~{}x)\leq E\mathopen{}\mathclose{{}\left(~{}x_{0}}\right)\}
\subset\mathbb{R}^{n}$.

Leaning on Proposition 3.7 and classical universal approximation results in [], it is further possible to establish the following error estimate also proven in Appendix A which gives an idea of the error accumulation rate that can be expected from this method.

###### Theorem 3.9 .

Suppose $~{}L,~{}M,E,S$ are nondegenerate metriplectic data such that $~{}L,~{}M$ have at least one continuous derivative, $E,S$ have Lipschitz continuous gradients, and at least one of $E,-S$ have bounded sublevel sets. For any $\varepsilon>0$, there exist nondegenerate metriplectic data $\tilde{~{}L},\tilde{~{}M},\tilde{E},\tilde{S}$ defined by two-layer neural networks such that, for all $T>0$,

$$ $\displaystyle\mathopen{}\mathclose{{}\left(\int_{0}^{T}\mathopen{}\mathclose{{ }\left|~{}x-\tilde{~{}x}}\right|^{2}\,dt}\right)^{\frac{1}{2}}\leq\varepsilon \mathopen{}\mathclose{{}\left|\frac{b}{a}}\right|\sqrt{e^{2aT}-2e^{aT}+T+1},$ $$

where $a,b\in\mathbb{R}$ are constants depending on both sets of metriplectic data and $\dot{\tilde{~{}x}}=\tilde{~{}L}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}
\right)\nabla\tilde{E}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)+
\tilde{~{}M}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)\nabla\tilde{S}
\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)$.

###### Remark 3.10 .

Theorem 3.9 provides a bound on state error over time under the assumption that the approximation error in the metriplectic networks can be controlled. On the other hand, notice that Theorem 3.9 can also be understood as a generic error bound on any trained metriplectic networks $\tilde{~{}L},\tilde{~{}M},\tilde{E},\tilde{S}$ provided universal approximation results are not invoked in the estimation leading to $\varepsilon b$.

This result confirms that the error in the state $~{}x$ for a fixed final time $T$ tends to zero with the approximation error in the networks $\tilde{~{}L},\tilde{~{}M},\tilde{E},\tilde{S}$, as one would hope based on the approximation capabilities of neural networks. More importantly, Theorem 3.9 also bounds the generalization error of any trained metriplectic network for an appropriate (and possibly large) $\varepsilon$ equal to the maximum approximation error on $K$, where the learned metriplectic trajectories are confined for all time. With this theoretical guidance, the remaining goal of this work is to demonstrate that NMS is also practically effective at learning metriplectic systems from data and exhibits reasonable generalization to unseen timescales.

###### Proposition 3.7 .

###### Remark 3.8 .

###### Theorem 3.9 .

###### Remark 3.10 .

## 4 Algorithm

Similar to previous approaches in [] and [], the learnable parameters in NMS are calibrated using data along solution trajectories to a given dynamical system. This brings up an important question regarding how much information about the system in question is realistically present in the training data. While many systems have a known metriplectic form, it is not always the case that one will know metriplectic governing equations for a given set of training data. Therefore, two approaches are considered in the experiments below corresponding to whether full or partial state information is assumed available during NMS training. More precisely, the state $~{}x=(~{}x^{o},~{}x^{u})$ will be partitioned into “observable” and “unobservable” variables, where $~{}x^{u}$ may be empty in the case that full state information is available. In a partially observable system $~{}x^{o}$ typically contains positions and momenta while $~{}x^{u}$ contains entropy or other configuration variables which are more difficult to observe during physical experiments. In both cases, NMS will learn a metriplectic system in $~{}x$ according to the procedure described in Algorithm 1.

Figure: Algorithm 1 Training neural metriplectic systems

Note that the batch-wise training strategy in Algorithm 1 requires initial conditions for $~{}x^{u}$ in the partially observed case. There are several options for this, and two specific strategies will be considered here. Suppose the data are drawn from the training interval $[0,T]$ with initial state $~{}x_{0}$ and final state $~{}x_{T}$. The first strategy sets $~{}x_{0}^{u}=~{}0$, $~{}x_{T}^{u}=~{}1$ (where $~{}1$ is the all ones vector), and $~{}x_{s}^{u}=~{}1/T$, $0<s<T$, so that the unobserved states are initially assumed to lie on a straight line. The second strategy is more sophisticated, and involves training a diffusion model to predict the distribution of $~{}x^{u}$ given $~{}x^{o}$. Specific details of this procedure are given in Appendix E.

## 5 Examples

The goal of the following experiments is to show that NMS is effective even when entropic information cannot be observed during training, yielding superior performance when compared to previous methods including GNODE, GFINN, and SPNN discussed in Section 2. The metrics considered for this purpose will be mean absolute error (MAE) and mean squared error (MSE) defined in the usual way as the average $\ell^{1}$ (resp. squared $\ell^{2}$) error between the discrete states $~{}x,\tilde{~{}x}\in\mathbb{R}^{n\times n_{s}}$. For brevity, many implementation details have been omitted here and can be found in Appendix B. An additional experiment showing the effectiveness of NMS in the presence of both full and partial state information can be found in Appendix C.

###### Remark 5.1 .

To facilitate a more equal parameter count between the compared metriplectic methods, the results of the experiments below were generated using the alternative parameterization $~{}D=~{}K~{}K^{\intercal}$ where $~{}K:K\to\mathbb{R}^{r\times r^{\prime}}$ is full and $r^{\prime}\geq r$. Of course, this change does not affect metriplecticity since $~{}D$ is still positive semi-definite for each $~{}x\in K$.

###### Remark 5.1 .

### 5.1 Two gas containers

The first test of NMS involves two ideal gas containers separated by movable wall which is removed at time $t_{0}$, allowing for the exchange of heat and volume. In this example, the motion of the state $~{}x=\begin{pmatrix}q&p&S_{1}&S_{2}\end{pmatrix}^{\intercal}$ is governed by the metriplectic equations:

$$ $\displaystyle\dot{q}$ $\displaystyle=\frac{p}{m},$ $\displaystyle\dot{p}=\frac{2}{3}\mathopen{}\mathclose{{}\left(\frac{E_{1}(~{}x )}{q}-\frac{E_{2}(~{}x)}{2L-q}}\right),$ $\displaystyle\dot{S}_{1}$ $\displaystyle=\frac{9N^{2}k_{B}^{2}\alpha}{4E_{1}(~{}x)}\mathopen{}\mathclose{ {}\left(\frac{1}{E_{1}(~{}x)}-\frac{1}{E_{2}(~{}x)}}\right),\quad$ $\displaystyle\dot{S}_{2}=-\frac{9N^{2}k_{B}^{2}\alpha}{4E_{1}(~{}x)}\mathopen{ }\mathclose{{}\left(\frac{1}{E_{1}(~{}x)}-\frac{1}{E_{2}(~{}x)}}\right),$ $$

where $(q,p)$ are the position and momentum of the separating wall, $S_{1},S_{2}$ are the entropies of the two subsystems, and the internal energies $E_{1},E_{2}$ are determined from the Sackur-Tetrode equation for ideal gases,
$S_{i}/Nk_{B}=\ln\mathopen{}\mathclose{{}\left(\hat{c}V_{i}E_{i}^{3/2}}\right),
1\leq i\leq 2.$
Here, $m$ denotes the mass of the wall, $2L$ is the total length of the system, and $V_{i}$ is the volume of the $i^{\mathrm{th}}$ container. As in [, ] $Nk_{B}=1$ and $\alpha=0.5$ fix the characteristic macroscopic unit of entropy while $\hat{c}=102.25$ ensures the argument of the logarithm defining $E_{i}$ is dimensionless. This leads to the total entropy $S(~{}x)=S_{1}+S_{2}$ and the total energy $E(~{}x)=(1/2m)p^{2}+E_{1}(~{}x)+E_{2}(~{}x),$ which are guaranteed to be nondecreasing and constant, respectively.

The primary goal here is to verify that NMS can accurately and stably predict gas container dynamics without the need to observe the entropic variables $S_{1},S_{2}$.
To that end, NMS has been compared to GNODE, SPNN, and GFINN on the task of predicting the trajectories of this metriplectic system over time, with results displayed in Table 1. More precisely, given an intial condition $~{}x_{0}$ and an interval $0<t_{\mathrm{train}}<t_{\mathrm{valid}}<t_{\mathrm{test}}$, each method is trained on partial state information (in the case of NMS) or full state information (in the case of the others) from the interval $[0,t_{\mathrm{train}}]$ and validated on $(t_{\mathrm{train}},t_{\mathrm{valid}}]$ before state errors in $q,p$ only are calculated on the whole interval $[0,t_{\mathrm{test}}]$. As can be seen from Table 1 and Figure 1, NMS is remarkably accurate over unseen timescales even in this unfair comparison, avoiding the unphysical behavior which often hinders soft-constrained methods like SPNN. The energy and instantaneous entropy plots in Figure 1 further confirm that the strong enforcement of metriplectic structure guaranteed by NMS leads to correct energetic and entropic dynamics for all time.

Figure: (a) Hamiltonian $H=\frac{p^{2}}{2m}$
Refer to caption: extracted/5627690/figs/tgc_legend.png

### 5.2 Thermoelastic double pendulum

Next, consider the thermoelastic double pendulum from [] with 10-dimensional state variable $~{}x=\begin{pmatrix}~{}q_{1}&~{}q_{2}&~{}p_{1}&~{}p_{2}&S_{1}&S_{2}\end{
pmatrix}^{\intercal}$, which represents a highly challenging benchmark for metriplectic methods. The equations of motion in this case are given for $1\leq i\leq 2$ as

$$ $\displaystyle\dot{~{}q}_{i}$ $\displaystyle=\frac{~{}p_{i}}{m_{i}},\quad\dot{~{}p}_{i}=-\partial_{~{}q_{i}} \mathopen{}\mathclose{{}\left(E_{1}(~{}x)+E_{2}(~{}x)}\right),\quad\dot{S}_{1} =\kappa\mathopen{}\mathclose{{}\left(T_{1}^{-1}T_{2}-1}\right),\quad\dot{S}_{2 }=\kappa\mathopen{}\mathclose{{}\left(T_{1}T_{2}^{-1}-1}\right),$ $$

where $\kappa>0$ is a thermal conductivity constant (set to 1), $m_{i}$ is the mass of the $i^{\mathrm{th}}$ spring (also set to 1) and $T_{i}=\partial_{S_{i}}E_{i}$ is its absolute temperature. In this case, $~{}q_{i},~{}p_{i}\in\mathbb{R}^{2}$ represent the position and momentum of the $i^{\mathrm{th}}$ mass, while $S_{i}$ represents the entropy of the $i^{\mathrm{th}}$ pendulum. As before, the total entropy $S(~{}x)=S_{1}+S_{2}$ is the sum of the entropies of the two springs, while defining the internal energies $E_{i}(~{}x)=(1/2)\mathopen{}\mathclose{{}\left(\ln\lambda_{i}}\right)^{2}+\ln
\lambda_{i}+e^{S_{i}-\ln\lambda_{i}}-1,\lambda_{1}=\mathopen{}\mathclose{{}
\left|~{}q_{i}}\right|,\lambda_{2}=\mathopen{}\mathclose{{}\left|~{}q_{2}-~{}q
_{1}}\right|,$
leads to the total energy $E(~{}x)=(1/2m_{1})\mathopen{}\mathclose{{}\left|~{}p_{1}}\right|^{2}+(1/2m_{2}
)\mathopen{}\mathclose{{}\left|~{}p_{2}}\right|^{2}+E_{1}(~{}x)+E_{2}(~{}x).$

**Table 1: Prediction errors for $~{}x^{o}$ measured in MSE and MAE on the interval $[0,t_{\mathrm{test}}]$ in the two gas containers example (left) and on the test set in the thermoelastic double pendulum example (right).**
|  | NODE | SPNN | GNODE | GFINN | NMS |
| --- | --- | --- | --- | --- | --- |
| MSE | .12 $\pm$ .04 | .13 $\pm$ .10 | .16 $\pm$ .10 | .07 $\pm$ .03 | .01 $\pm$ .02 |
| MAE | .25 $\pm$ .10 | .26 $\pm$ .14 | .25 $\pm$ .13 | .13 $\pm$.03 | .08 $\pm$ .06 |

The task in this case is prediction across initial conditions. As in [], 100 trajectories are drawn from the ranges in Appendix B and integrated over the interval $[0,40]$ with $\Delta t=0.1$, with an 80/10/10 split for training/validation/testing. Here all compared models are trained using full state information. As seen in Table 1, NMS is again the most performant, although all models struggle to approximate the dynamics over the entire training interval. It is also notable that the training time of NMS is greatly decreased relative to GNODE and GFINN due to its improved quadratic scaling; a representative study to this effect is given in Appendix D.

## 6 Conclusion

Neural metriplectic systems (NMS) have been considered for learning finite-dimensional metriplectic dynamics from data. Making use of novel non-redundant parameterizations for metriplectic operators, NMS provably approximates arbitrary nondegenerate metriplectic systems with generalization error bounded in terms of the operator approximation quality. Benchmark examples have shown that NMS is both more scalable and more accurate than previous methods, including when only partial state information is observed. Future work will consider extensions of NMS to infinite-dimensional metriplectic systems with the aim of addressing its main limitation: the difficulty of scaling NMS (among all present methods for metriplectic learning) to realistic, 3-D problems of the size that would be considered in practice. A promising direction is to consider the use of NMS in model reduction, where sparse, large-scale systems are converted to small, dense systems through a clever choice of encoding/decoding.

## Acknowledgements

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. This article has been co-authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns all right, title and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-plan. The work of A. Gruber and N. Trask is supported by the John von Neumann fellowship at Sandia, the U.S. Department of Energy, Office of Advanced Computing Research under the “Scalable and Efficient Algorithms - Causal Reasoning, Operators, Graphs and Spikes" (SEA-CROGS) project, and the DoE Early Career Research Program. K. Lee acknowledges support from the U.S. National Science Foundation under grant IIS 2338909. N. Park and H. Lim acknowledge support from the Korea Advanced Institute of Science and Technology (KAIST) grant funded by the Korea government (MSIT) (No. G04240001, Physics-inspired Deep Learning), as well as the Institute for Information & Communications Technology Planning & Evaluation (IITP) grants funded by the Korea government (MSIT) (No. 2020-0-01361, Artificial Intelligence Graduate School Program (Yonsei University)). The authors thank Jonas Actor for careful proofreading and suggestions regarding the presentation of Appendix A.

## Appendix A Proof of Theoretical Results

This Appendix provides proof of the analytical results in Section 3 of the body. First, the parameterizations of $~{}L,~{}M$ in terms of exterior algebra are established.

###### Proof of Lemma 3.2 .

First, it is necessary to check that the operators $~{}L,~{}M$ parameterized this way satisfy the symmetries and degeneracy conditions claimed in the statement. To that end, recall that $~{}a\wedge~{}b\simeq~{}a~{}b^{\intercal}-~{}b~{}a^{\intercal}$, meaning that $\mathopen{}\mathclose{{}\left(~{}a~{}b^{\intercal}-~{}b~{}a^{\intercal}}\right
)^{\intercal}\simeq~{}b\wedge~{}a=-~{}a\wedge~{}b$. It follows that $~{}A^{\intercal}\simeq\tilde{\mathsf{A}}=-\mathsf{A}$ where $\tilde{\mathsf{A}}$ denotes the reversion of $\mathsf{A}$, i.e., $\tilde{\mathsf{A}}=\sum_{i<j}A^{ij}~{}e_{j}\wedge~{}e_{i}.$ Therefore, we may write

$$ $\displaystyle~{}L^{\intercal}$ $\displaystyle\simeq\tilde{\mathsf{A}}-\frac{1}{\mathopen{}\mathclose{{}\left| \nabla S}\right|^{2}}\ThisStyle{\stackengine{-.1\LMpt}{$\SavedStyle~{}A\nabla S \wedge\nabla S$}{\stretchto{\scaleto{\SavedStyle\mkern 0.2mu\AC}{.5150}}{.6}}{ O}{c}{F}{T}{S}}=-\mathsf{A}+\frac{1}{\mathopen{}\mathclose{{}\left|\nabla S} \right|^{2}}~{}A\nabla S\wedge\nabla S\simeq-~{}L$ $$

showing that $~{}L^{\intercal}=-~{}L$. Moreover, using that

$$ $\mathopen{}\mathclose{{}\left(~{}b\wedge~{}c}\right)\cdot~{}a=-~{}a\cdot \mathopen{}\mathclose{{}\left(~{}b\wedge~{}c}\right)=\mathopen{}\mathclose{{} \left(~{}a\cdot~{}c}\right)~{}b-\mathopen{}\mathclose{{}\left(~{}a\cdot~{}b} \right)~{}c,$ $$

it follows that

$$ $\displaystyle~{}L\nabla S$ $\displaystyle=\mathsf{A}\cdot\nabla S-\frac{1}{\mathopen{}\mathclose{{}\left| \nabla S}\right|^{2}}\mathopen{}\mathclose{{}\left(~{}A\nabla S\wedge\nabla S} \right)\cdot\nabla S=~{}A\nabla S-~{}A\nabla S=~{}0,$ $$

since $\nabla S\cdot~{}A\nabla S=-\nabla S\cdot~{}A\nabla S=0$. Moving to the case of $~{}M$, notice that $~{}M=D_{st}~{}v^{s}\otimes~{}v^{t}$ for a particular choice of $~{}v$, meaning that

$$ $\displaystyle~{}M^{\intercal}$ $\displaystyle=\sum_{s,t}D_{st}\mathopen{}\mathclose{{}\left(~{}v^{s}\otimes~{} v^{t}}\right)^{\intercal}=\sum_{s,t}D_{st}~{}v^{t}\otimes~{}v^{s}=\sum_{t,s}D_ {ts}~{}v^{s}\otimes~{}v^{t}=\sum_{s,t}D_{st}~{}v^{s}\otimes~{}v^{t}=~{}M,$ $$

since $~{}D$ is a symmetric matrix. Additionally, it is straightforward to check that, for any $1\leq s\leq r$,

$$ $\displaystyle~{}v^{s}\cdot\nabla E$ $\displaystyle=\mathopen{}\mathclose{{}\left(~{}b^{s}-\frac{~{}b^{s}\cdot\nabla E }{\mathopen{}\mathclose{{}\left|\nabla E}\right|^{2}}\nabla E}\right)\cdot \nabla E=~{}b^{s}\cdot\nabla E-~{}b^{s}\cdot\nabla E=0.$ $$

So, it follows immediately that

$$ $\displaystyle~{}M\nabla E$ $\displaystyle=\sum_{s,t}D_{st}\mathopen{}\mathclose{{}\left(~{}v^{s}\otimes~{} v^{t}}\right)\cdot\nabla E=\sum_{s,t}D_{st}\mathopen{}\mathclose{{}\left(~{}v^ {t}\cdot\nabla E}\right)~{}v^{s}=~{}0.$ $$

Now, observe that

$$ $\displaystyle~{}L$ $\displaystyle=~{}A-\frac{1}{\mathopen{}\mathclose{{}\left|\nabla S}\right|^{2} }\mathopen{}\mathclose{{}\left(~{}A\nabla S\mathopen{}\mathclose{{}\left( \nabla S}\right)^{\intercal}-\nabla S\mathopen{}\mathclose{{}\left(~{}A\nabla S }\right)^{\intercal}}\right)$ $\displaystyle=~{}A-\frac{1}{\mathopen{}\mathclose{{}\left|\nabla S}\right|^{2} }\mathopen{}\mathclose{{}\left(~{}A\nabla S\mathopen{}\mathclose{{}\left( \nabla S^{\intercal}}\right)+\nabla S\mathopen{}\mathclose{{}\left(\nabla S} \right)^{\intercal}~{}A}\right)$ $\displaystyle=\mathopen{}\mathclose{{}\left(~{}I-\frac{\nabla S\mathopen{} \mathclose{{}\left(\nabla S}\right)^{\intercal}}{\mathopen{}\mathclose{{}\left |\nabla S}\right|^{2}}}\right)~{}A\mathopen{}\mathclose{{}\left(~{}I-\frac{ \nabla S\mathopen{}\mathclose{{}\left(\nabla S}\right)^{\intercal}}{\mathopen{ }\mathclose{{}\left|\nabla S}\right|^{2}}}\right)=~{}P^{\perp}_{S}~{}A~{}P^{ \perp}_{S},$ $$

since $~{}A^{\intercal}=-~{}A$ and hence $~{}v^{\intercal}~{}A~{}v=0$ for all $~{}v\in\mathbb{R}^{n}$. Similarly, it follows that for every $1\leq s\leq r$,

$$ $\displaystyle~{}P^{\perp}_{E}~{}b^{s}=~{}b^{s}-\frac{~{}b^{s}\cdot\nabla E}{ \mathopen{}\mathclose{{}\left|\nabla E}\right|^{2}}\nabla E,$ $$

and therefore $~{}M$ is expressible as

$$ $\displaystyle~{}M=\sum_{s,t}D_{st}\mathopen{}\mathclose{{}\left(~{}P^{\perp}_{ E}~{}b^{s}}\right)\mathopen{}\mathclose{{}\left(~{}P^{\perp}_{E}~{}b^{t}} \right)^{\intercal}=~{}P^{\perp}_{E}~{}B~{}D~{}B^{\intercal}~{}P^{\perp}_{E}.\quad\qed$ $$

With Lemma 3.2 established, the proof of Theorem 3.4 is straightforward.

###### Proof of Theorem 3.4 .

The “if” direction follows immediately from Lemma 3.2. Now, suppose that $~{}L$ and $~{}M$ define a metriplectic system, meaning that the mentioned symmetries and degeneracy conditions hold. Then, it follows from $~{}L\nabla S=~{}0$ that the projection $~{}P^{\perp}_{S}~{}L~{}P^{\perp}_{S}=~{}L$ leaves $~{}L$ invariant, so that choosing $~{}A=~{}L$ yields $~{}P^{\perp}_{S}~{}A~{}P^{\perp}_{S}=~{}L$. Similarly, from positive semi-definiteness and $~{}M\nabla E=~{}0$ it follows that $~{}M=~{}U~{}\Lambda~{}U^{\intercal}=~{}P^{\perp}_{E}~{}U~{}\Lambda~{}U^{
\intercal}~{}P^{\perp}_{E}$ for some column-orthonormal $~{}U\in\mathbb{R}^{N\times r}$ and positive diagonal $~{}\Lambda\in\mathbb{R}^{r\times r}$. Therefore, choosing $~{}B=~{}U$ and $~{}D=~{}\Lambda$ yields $~{}M=~{}P^{\perp}_{E}~{}B~{}D~{}B^{\intercal}~{}P^{\perp}$, as desired.
∎

Looking toward the proof of Proposition 3.7, we also need to establish the following Lemmata which give control over the orthogonal projectors $~{}P_{\tilde{E}}^{\perp},~{}P_{\tilde{S}}^{\perp}$. First, we recall how control over the $L^{\infty}$ norm $\mathopen{}\mathclose{{}\left|\cdot}\right|_{\infty}$ of a matrix field gives control over its spectral norm $\mathopen{}\mathclose{{}\left|\cdot}\right|$.

###### Lemma A.1 .

Let $~{}A:K\to\mathbb{R}^{n\times n}$ be a matrix field defined on the compact set $K\subset\mathbb{R}^{n}$ with $m$ continuous derivatives. Then, for any $\varepsilon>0$ there exists a two-layer neural network $\tilde{~{}A}:K\to\mathbb{R}^{n\times n}$ such that $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}A-\tilde{~{}A}}\right|<\varepsilon$ and $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|\nabla^{k}~{}A-\nabla^{k}\tilde{
~{}A}}\right|_{\infty}<\varepsilon$ for $1\leq k\leq m$ where $\nabla^{k}$ is the (total) derivative operator of order $k$.

###### Proof.

This will be a direct consequence of Corollary 2.2 in [] provided we show that $\mathopen{}\mathclose{{}\left|~{}A}\right|\leq c\mathopen{}\mathclose{{}\left|
~{}A}\right|_{\infty}$ for some $c>0$. To that end, if $\sigma_{1}\geq...\geq\sigma_{r}>0$ ($r\leq n$) denote the nonzero singular values of $~{}A-\tilde{~{}A}$, it follows that for each $~{}x\in K$,

$$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}A-\tilde{~{}A}}\right|$ $\displaystyle=\sigma_{1}\leq\sqrt{\sigma_{1}^{2}+...+\sigma_{r}^{2}}=\sqrt{ \sum_{i,j}\mathopen{}\mathclose{{}\left|A_{ij}-\tilde{A}_{ij}}\right|^{2}}= \mathopen{}\mathclose{{}\left|~{}A-\tilde{~{}A}}\right|_{F}.$ $$

On the other hand, it also follows that

$$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}A-\tilde{~{}A}}\right|_{F}$ $\displaystyle=\sqrt{\sum_{i,j}\mathopen{}\mathclose{{}\left|A_{ij}-\tilde{A}_{ ij}}\right|^{2}}\leq\sqrt{\sum_{i,j}\max_{i,j}\mathopen{}\mathclose{{}\left|A_ {ij}-\tilde{A}_{ij}}\right|}=n\sqrt{\max_{i,j}\mathopen{}\mathclose{{}\left|A_ {ij}-\tilde{A}_{ij}}\right|}=n\mathopen{}\mathclose{{}\left|~{}A-\tilde{~{}A}} \right|_{\infty},$ $$

and therefore the desired inequality holds with $c=n$. Now, for any $\varepsilon>0$ it follows from [] that there exists a two layer network $\tilde{~{}A}$ with $m$ continuous derivatives such that $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}A-\tilde{~{}A}}\right|_{
\infty}<\varepsilon/n$ and $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|\nabla^{k}~{}A-\nabla^{k}\tilde{
~{}A}}\right|_{\infty}<\varepsilon/n<\varepsilon$ for all $1\leq k\leq m$. Therefore, it follows that

$$ $\displaystyle\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}A-\tilde{~{}A}} \right|\leq n\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}A-\tilde{~{}A}} \right|_{\infty}<n\frac{\varepsilon}{n}=\varepsilon,$ $$

completing the argument.
∎

Next, we bound the deviation in the orthogonal projectors $~{}P_{\tilde{E}}^{\perp},~{}P_{\tilde{S}}^{\perp}$.

###### Lemma A.2 .

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be such that $\nabla f\neq~{}0$ on the compact set $K\subset\mathbb{R}^{n}$. For any $\varepsilon>0$, there exists a two-layer neural network $\tilde{f}:K\to\mathbb{R}$ such that $\nabla\tilde{f}\neq~{}0$ on $K$, $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|f-\tilde{f}}\right|<\varepsilon,
\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|\nabla f-\nabla\tilde{f}}\right|<\varepsilon$, and $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}P_{f}^{\perp}-~{}P_{\tilde{f}
}^{\perp}}\right|<\varepsilon$.

###### Proof.

Denote $\nabla f=~{}v$ and consider any $\tilde{~{}v}:K\to\mathbb{R}$. Since $\mathopen{}\mathclose{{}\left|~{}v}\right|\leq\mathopen{}\mathclose{{}\left|
\tilde{~{}v}}\right|+\mathopen{}\mathclose{{}\left|~{}v-\tilde{~{}v}}\right|$, it follows for all $~{}x\in K$ that whenever $\mathopen{}\mathclose{{}\left|~{}v-\tilde{~{}v}}\right|<(1/2)\inf_{~{}x\in K}
\mathopen{}\mathclose{{}\left|~{}v}\right|$,

$$ $\mathopen{}\mathclose{{}\left|\tilde{~{}v}}\right|\geq\mathopen{}\mathclose{{} \left|~{}v}\right|-\mathopen{}\mathclose{{}\left|~{}v-\tilde{~{}v}}\right|> \mathopen{}\mathclose{{}\left|~{}v}\right|-\frac{1}{2}\inf_{~{}x\in K} \mathopen{}\mathclose{{}\left|~{}v}\right|>0,$ $$

so that $\tilde{~{}v}\neq 0$ in $K$, and since the square function is monotonic,

$$ $\inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|\tilde{~{}v}}\right|^{2}\geq\inf _{~{}x\in K}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left|~{}v} \right|-\frac{1}{2}\inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}v}\right|} \right)^{2}=\frac{1}{4}\inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}v} \right|^{2}.$ $$

On the other hand, we also have $\mathopen{}\mathclose{{}\left|\tilde{~{}v}}\right|\leq\mathopen{}\mathclose{{}
\left|~{}v}\right|+\mathopen{}\mathclose{{}\left|\tilde{~{}v}-~{}v}\right|<
\mathopen{}\mathclose{{}\left|~{}v}\right|+(1/2)\inf_{~{}x\in K}\mathopen{}
\mathclose{{}\left|~{}v}\right|$, so that, adding and subtracting $\tilde{~{}v}~{}v^{\intercal}$ and applying Cauchy-Schwarz, it follows that for all $~{}x\in K$,

$$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}v~{}v^{\intercal}-\tilde{~{}v} \tilde{~{}v}^{\intercal}}\right|$ $\displaystyle\leq\mathopen{}\mathclose{{}\left|~{}v-\tilde{~{}v}}\right| \mathopen{}\mathclose{{}\left|~{}v}\right|+\mathopen{}\mathclose{{}\left| \tilde{~{}v}}\right|\mathopen{}\mathclose{{}\left|~{}v-\tilde{~{}v}}\right| \leq 2\max\mathopen{}\mathclose{{}\left\{\mathopen{}\mathclose{{}\left|~{}v} \right|,\mathopen{}\mathclose{{}\left|\tilde{~{}v}}\right|}\right\}\mathopen{} \mathclose{{}\left|~{}v-\tilde{~{}v}}\right|<\mathopen{}\mathclose{{}\left(2 \mathopen{}\mathclose{{}\left|~{}v}\right|+\inf_{~{}x\in K}\mathopen{} \mathclose{{}\left|~{}v}\right|}\right)\mathopen{}\mathclose{{}\left|~{}v- \tilde{~{}v}}\right|.$ $$

Now, by Corollary 2.2 in [], for any $\varepsilon>0$ there exists a two-layer neural network $\tilde{f}:K\to\mathbb{R}$ such that

$$ $\displaystyle\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}v-\nabla\tilde{f }}\right|<\min\mathopen{}\mathclose{{}\left\{\frac{1}{2}\inf_{~{}x\in K} \mathopen{}\mathclose{{}\left|~{}v}\right|,\frac{\inf_{~{}x\in K}\mathopen{} \mathclose{{}\left|~{}v}\right|^{2}}{2\sup_{~{}x\in K}\mathopen{}\mathclose{{} \left|~{}v}\right|+\inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}v}\right|} \frac{\varepsilon}{4},\varepsilon}\right\}\leq\varepsilon,$ $$

and also $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|f-\tilde{f}}\right|<\varepsilon.$ Letting $\tilde{~{}v}=\nabla\tilde{f}$, it follows that for all $~{}x\in K$,

$$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}P_{f}^{\perp}-~{}P_{\tilde{f}}^{ \perp}}\right|$ $\displaystyle=\mathopen{}\mathclose{{}\left|\frac{~{}v~{}v^{\intercal}}{ \mathopen{}\mathclose{{}\left|~{}v}\right|^{2}}-\frac{\tilde{~{}v}\tilde{~{}v} ^{\intercal}}{\mathopen{}\mathclose{{}\left|\tilde{~{}v}}\right|^{2}}}\right| \leq\frac{\mathopen{}\mathclose{{}\left|~{}v~{}v^{\intercal}-\tilde{~{}v} \tilde{~{}v}^{\intercal}}\right|}{\min\mathopen{}\mathclose{{}\left\{\mathopen {}\mathclose{{}\left|~{}v}\right|^{2},\mathopen{}\mathclose{{}\left|\tilde{~{} v}}\right|^{2}}\right\}}\leq\frac{2\mathopen{}\mathclose{{}\left|~{}v}\right|+ \inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}v}\right|}{\min\mathopen{} \mathclose{{}\left\{\mathopen{}\mathclose{{}\left|~{}v}\right|^{2},\mathopen{} \mathclose{{}\left|\tilde{~{}v}}\right|^{2}}\right\}}\mathopen{}\mathclose{{} \left|~{}v-\tilde{~{}v}}\right|,$ $$

and therefore, taking the supremum of both sides and applying the previous work yields the desired estimate,

$$ $\displaystyle\sup_{~{}x\in K}$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}P_{f}^{\perp}-~{}P_{\tilde{f}}^{ \perp}}\right|\leq 4\frac{2\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}v} \right|+\inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}v}\right|}{\inf_{~{}x \in K}\mathopen{}\mathclose{{}\left|~{}v}\right|^{2}}\sup_{~{}x\in K}\mathopen {}\mathclose{{}\left|~{}v-\tilde{~{}v}}\right|<\varepsilon.\qed$ $$

With these intermediate results established, the proof of the approximation result Proposition 3.7 proceeds as follows.

###### Proof of Proposition 3.7 .

Recall from Theorem 3.4 that we can write $~{}L=~{}P^{\perp}_{S}\mathopen{}\mathclose{{}\left(~{}A_{\mathrm{tri}}-~{}A_{
\mathrm{tri}}^{\intercal}}\right)~{}P^{\perp}_{S}$ and similarly for $\tilde{~{}L}$. Notice that, by adding and subtracting $~{}P_{\tilde{S}}^{\perp}~{}A_{\mathrm{tri}}~{}P_{S}^{\perp}$ and $~{}P_{\tilde{S}}^{\perp}\tilde{~{}A}_{\mathrm{tri}}~{}P_{S}^{\perp}$, it follows that for all $~{}x\in K$,

$$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}P_{S}^{\perp}~{}A_{\mathrm{tri}} ~{}P_{S}^{\perp}-~{}P_{\tilde{S}}^{\perp}\tilde{~{}A}_{\mathrm{tri}}~{}P_{ \tilde{S}}^{\perp}}\right|$ $\displaystyle\qquad=\bigg{|}\mathopen{}\mathclose{{}\left(~{}P^{\perp}_{S}-~{} P^{\perp}_{\tilde{S}}}\right)~{}A_{\mathrm{tri}}~{}P^{\perp}_{S}+~{}P^{\perp}_ {\tilde{S}}\mathopen{}\mathclose{{}\left(~{}A_{\mathrm{tri}}-\tilde{~{}A}_{ \mathrm{tri}}}\right)~{}P^{\perp}_{S}+~{}P^{\perp}_{\tilde{S}}\tilde{~{}A}_{ \mathrm{tri}}\mathopen{}\mathclose{{}\left(~{}P^{\perp}_{S}-~{}P^{\perp}_{ \tilde{S}}}\right)\bigg{|}$ $\displaystyle\qquad\leq\mathopen{}\mathclose{{}\left|~{}P^{\perp}_{S}-~{}P^{ \perp}_{\tilde{S}}}\right|\mathopen{}\mathclose{{}\left|~{}A_{\mathrm{tri}}} \right|+\mathopen{}\mathclose{{}\left|~{}A_{\mathrm{tri}}-\tilde{~{}A}_{ \mathrm{tri}}}\right|+\mathopen{}\mathclose{{}\left|\tilde{~{}A}_{\mathrm{tri} }}\right|\mathopen{}\mathclose{{}\left|~{}P^{\perp}_{S}-~{}P^{\perp}_{\tilde{S }}}\right|$ $\displaystyle\qquad\leq 2\max\mathopen{}\mathclose{{}\left\{\mathopen{} \mathclose{{}\left|~{}A_{\mathrm{tri}}}\right|,\mathopen{}\mathclose{{}\left| \tilde{~{}A}_{\mathrm{tri}}}\right|}\right\}\mathopen{}\mathclose{{}\left|~{}P ^{\perp}_{S}-~{}P^{\perp}_{\tilde{S}}}\right|+\mathopen{}\mathclose{{}\left|~{ }A_{\mathrm{tri}}-\tilde{~{}A}_{\mathrm{tri}}}\right|$ $$

where we have used that $~{}P_{S}^{\perp},~{}P_{\tilde{S}}^{\perp}$ have unit spectral norm. By Lemma A.1, for any $\varepsilon>0$ there exists a two layer neural network $\tilde{~{}A}_{\mathrm{tri}}$ such that $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}A_{\mathrm{tri}}-\tilde{~{}A}
_{\mathrm{tri}}}\right|<\frac{\varepsilon}{4}$, and by Lemma A.2 there exists a two-layer network $\tilde{S}$ with $\nabla\tilde{S}\neq~{}0$ on $K$ such that

$$ $\displaystyle\sup_{~{}x\in K}$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}P_{S}^{\perp}-~{}P_{\tilde{S}}^{ \perp}}\right|<\min\mathopen{}\mathclose{{}\left\{\varepsilon,\max\mathopen{} \mathclose{{}\left\{\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}A_{ \mathrm{tri}}}\right|,\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|\tilde{~{} A}_{\mathrm{tri}}}\right|}\right\}^{-1}\frac{\varepsilon}{8}}\right\}.$ $$

It follows that $\tilde{S},\nabla\tilde{S}$ are $\varepsilon$-close to $S,\nabla S$ on $K$ and

$$ $\displaystyle\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left(2\max\mathopen{} \mathclose{{}\left\{\mathopen{}\mathclose{{}\left|~{}A_{\mathrm{tri}}}\right|, \mathopen{}\mathclose{{}\left|\tilde{~{}A}_{\mathrm{tri}}}\right|}\right\} \mathopen{}\mathclose{{}\left|~{}P^{\perp}_{S}-~{}P^{\perp}_{\tilde{S}}}\right |}\right)<\frac{\varepsilon}{4}.$ $$

Therefore, the estimate

$$ $\displaystyle\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}L-\tilde{~{}L}}\right|$ $\displaystyle\leq 2\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}P_{S}^{ \perp}~{}A_{\mathrm{tri}}~{}P_{S}^{\perp}-~{}P_{\tilde{S}}^{\perp}\tilde{~{}A} _{\mathrm{tri}}~{}P_{\tilde{S}}^{\perp}}\right|<2\mathopen{}\mathclose{{}\left (\frac{\varepsilon}{4}+\frac{\varepsilon}{4}}\right)=\varepsilon,$ $$

implies that $\tilde{~{}L}$ is $\varepsilon$-close to $~{}L$ on $K$ as well.

Moving to the case of $~{}M$, we see that for all $~{}x\in K$, by writing $~{}M=~{}U~{}\Lambda~{}U^{\intercal}=~{}K_{\mathrm{chol}}~{}K_{\mathrm{chol}}^{\intercal}$ for $~{}K_{\mathrm{chol}}=~{}U~{}\Lambda^{1/2}$ and repeating the first calculation with $~{}K_{\mathrm{chol}}$ in place of $~{}A_{\mathrm{tri}}$ and $~{}P_{E}^{\perp}$ in place of $~{}P_{S}^{\perp}$,

$$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}P_{E}^{\perp}~{}K_{\mathrm{chol} }~{}K_{\mathrm{chol}}^{\intercal}~{}P_{E}^{\perp}-~{}P_{\tilde{E}}^{\perp} \tilde{~{}K}_{\mathrm{chol}}\tilde{~{}K}_{\mathrm{chol}}^{\intercal}~{}P_{ \tilde{E}}^{\perp}}\right|$ $\displaystyle\qquad\leq 2\max\mathopen{}\mathclose{{}\left\{\mathopen{} \mathclose{{}\left|~{}K_{\mathrm{chol}}}\right|,\mathopen{}\mathclose{{}\left| \tilde{~{}K}_{\mathrm{chol}}}\right|}\right\}\mathopen{}\mathclose{{}\left|~{} P^{\perp}_{E}-~{}P^{\perp}_{\tilde{E}}}\right|+\mathopen{}\mathclose{{}\left|~ {}K_{\mathrm{chol}}~{}K_{\mathrm{chol}}^{\intercal}-\tilde{~{}K}_{\mathrm{chol }}\tilde{~{}K}_{\mathrm{chol}}^{\intercal}}\right|.$ $$

Moreover, if $\mathopen{}\mathclose{{}\left|~{}K_{\mathrm{chol}}-\tilde{~{}K}_{\mathrm{chol}
}}\right|<(1/2)\inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}K_{\mathrm{
chol}}}\right|$ for all $~{}x\in K$ then similar arguments as used in the proof of Lemma A.2 yield the following estimate for all $~{}x\in K$,

$$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}K_{\mathrm{chol}}~{}K_{\mathrm{ chol}}^{\intercal}-\tilde{~{}K}_{\mathrm{chol}}\tilde{~{}K}_{\mathrm{chol}}^{ \intercal}}\right|$ $\displaystyle\leq 2\max\mathopen{}\mathclose{{}\left\{\mathopen{}\mathclose{{} \left|~{}K_{\mathrm{chol}}}\right|,\mathopen{}\mathclose{{}\left|\tilde{~{}K}_ {\mathrm{chol}}}\right|}\right\}\mathopen{}\mathclose{{}\left|~{}K_{\mathrm{ chol}}-\tilde{~{}K}_{\mathrm{chol}}}\right|$ $\displaystyle\leq\mathopen{}\mathclose{{}\left(2\mathopen{}\mathclose{{}\left| ~{}K_{\mathrm{chol}}}\right|+\inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{} K_{\mathrm{chol}}}\right|}\right)\mathopen{}\mathclose{{}\left|~{}K_{\mathrm{ chol}}-\tilde{~{}K}_{\mathrm{chol}}}\right|.$ $$

As before, we now invoke Lemma A.1 to construct a two-layer lower-triangular network $\tilde{~{}K}_{\mathrm{chol}}$ such that

$$ $\displaystyle\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}K_{\mathrm{chol} }-\tilde{~{}K}_{\mathrm{chol}}}\right|<\min\bigg{\{}\frac{1}{2}\inf_{~{}x\in K }\mathopen{}\mathclose{{}\left|~{}K_{\mathrm{chol}}}\right|,\mathopen{} \mathclose{{}\left(2\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}K_{ \mathrm{chol}}}\right|+\inf_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}K_{ \mathrm{chol}}}\right|}\right)^{-1}\frac{\varepsilon}{2}\bigg{\}},$ $$

as well as (using Lemma A.2) a network $\tilde{E}$ satisfying $\nabla\tilde{E}\neq~{}0$ on $K$ and

$$ $\displaystyle\sup_{~{}x\in K}$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}P_{E}^{\perp}-~{}P_{\tilde{E}}^{ \perp}}\right|<\min\mathopen{}\mathclose{{}\left\{\varepsilon,\max\mathopen{} \mathclose{{}\left\{\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}K_{ \mathrm{chol}}}\right|,\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|\tilde{~{ }K}_{\mathrm{chol}}}\right|}\right\}^{-1}\frac{\varepsilon}{4}}\right\}.$ $$

Again, it follows that $\tilde{E},\nabla\tilde{E}$ are $\varepsilon$-close to $E,\nabla E$ on $K$, and by the work above we conclude

$$ $\displaystyle\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}M-\tilde{~{}M}}\right|$ $\displaystyle=\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}P_{E}^{\perp}~{ }K_{\mathrm{chol}}~{}K_{\mathrm{chol}}^{\intercal}~{}P_{E}^{\perp}-~{}P_{ \tilde{E}}^{\perp}\tilde{~{}K}_{\mathrm{chol}}\tilde{~{}K}_{\mathrm{chol}}^{ \intercal}~{}P_{\tilde{E}}^{\perp}}\right|<\frac{\varepsilon}{2}+\frac{ \varepsilon}{2}=\varepsilon,$ $$

as desired.
∎

It is now possible to give a proof of the error bound in Theorem 3.9. Recall the $L^{2}\mathopen{}\mathclose{{}\left([0,T]}\right)$ error metric $\mathopen{}\mathclose{{}\left\|~{}x}\right\|$ and Lipschitz constant $L_{f}$, defined for all $~{}x,~{}y\in\mathbb{R}^{n}$ and Lipschitz continuous functions $f$ as

$$ $\displaystyle\mathopen{}\mathclose{{}\left\|~{}x}\right\|^{2}=\int_{0}^{T} \mathopen{}\mathclose{{}\left|~{}x}\right|^{2}\,dt,\quad\mathopen{}\mathclose{ {}\left|f(~{}x)-f(~{}y)}\right|\leq L_{f}\mathopen{}\mathclose{{}\left|~{}x-~{ }y}\right|.$ $$

###### Proof of Theorem 3.9 .

First, note that the assumption that one of $E,-S$ (without loss of generality, say $E$) has bounded sublevel sets implies bounded trajectories for the state $~{}x$ as in Remark 3.8, so we may assume $~{}x\in K$ for some compact $K\subset\mathbb{R}^{n}$. Moreover, for any $\varepsilon>0$ it follows from Proposition 3.7 that there are approximate networks $\tilde{E},\tilde{S}$ which are $\varepsilon$-close to $E,S$ on $K$. Additionally, it follows that $\tilde{E},\tilde{S}$ have nonzero gradients $\nabla\tilde{E},\nabla\tilde{S}$ which are also $\varepsilon$-close to the true gradients $\nabla E,\nabla S$ on $K$. This implies that for each $~{}x\in K$, $E=\tilde{E}+(E-\tilde{E})\leq\tilde{E}+\varepsilon$, so it follows that the sublevel sets $\{~{}x\,|\,\tilde{E}(~{}x)\leq m\}\subseteq\{~{}x\,|\,E(~{}x)\leq m+\varepsilon\}$ are also bounded. Therefore, we may assume (by potentially enlarging $K$) that both $~{}x,\tilde{~{}x}\in K$ lie in the compact set $K$ for all time.

Now, let $~{}y=~{}x-\tilde{~{}x}$. The next goal is to bound the following quantity:

$$ $\displaystyle\mathopen{}\mathclose{{}\left|\dot{~{}y}}\right|$ $\displaystyle=\bigg{|}~{}L(~{}x)\nabla E(~{}x)+~{}M(~{}x)\nabla S(~{}x)-\tilde {~{}L}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)\nabla\tilde{E} \mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)-\tilde{~{}M}\mathopen{} \mathclose{{}\left(\tilde{~{}x}}\right)\nabla\tilde{S}\mathopen{}\mathclose{{} \left(\tilde{~{}x}}\right)\bigg{|}$ $\displaystyle=\bigg{|}\mathopen{}\mathclose{{}\left(~{}L(~{}x)\nabla E(~{}x)- \tilde{~{}L}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)\nabla E \mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)}\right)+\mathopen{} \mathclose{{}\left(~{}M(~{}x)\nabla S(~{}x)-\tilde{~{}M}\mathopen{}\mathclose{ {}\left(\tilde{~{}x}}\right)\nabla S\mathopen{}\mathclose{{}\left(\tilde{~{}x} }\right)}\right)\bigg{|}=:\mathopen{}\mathclose{{}\left|\dot{~{}y}_{E}+\dot{~{ }y}_{S}}\right|.$ $$

To that end, notice that by adding and subtracting $~{}L(~{}x)\nabla E\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right),\tilde{~{
}L}(~{}x)\nabla E\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right),\tilde{~{}
L}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)\nabla E\mathopen{}
\mathclose{{}\left(\tilde{~{}x}}\right)$, it follows that

$$ $\displaystyle\dot{~{}y}_{E}$ $\displaystyle=~{}L(~{}x)\mathopen{}\mathclose{{}\left(\nabla E(~{}x)-\nabla E \mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)}\right)+\mathopen{} \mathclose{{}\left(~{}L(~{}x)-\tilde{~{}L}(~{}x)}\right)\nabla E\mathopen{} \mathclose{{}\left(\tilde{~{}x}}\right)$ $\displaystyle\qquad+\mathopen{}\mathclose{{}\left(\tilde{~{}L}(~{}x)-\tilde{~{ }L}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)}\right)\nabla E\mathopen {}\mathclose{{}\left(\tilde{~{}x}}\right)+\tilde{~{}L}\mathopen{}\mathclose{{} \left(\tilde{~{}x}}\right)\mathopen{}\mathclose{{}\left(\nabla E\mathopen{} \mathclose{{}\left(\tilde{~{}x}}\right)-\nabla\tilde{E}\mathopen{}\mathclose{{ }\left(\tilde{~{}x}}\right)}\right).$ $$

By Proposition 3.7 there exists a two-layer neural network $\tilde{~{}L}$ with one continuous derivative such that $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}L-\tilde{~{}L}}\right|<\varepsilon$, which implies that $\tilde{~{}L}$ is Lipschitz continuous with (uniformly well-approximated) Lipschitz constant. Using this fact along with the assumed Lipschitz continuity of $\nabla E$ and the approximation properties of the network $\tilde{E}$ already constructed then yields

$$ $\displaystyle\mathopen{}\mathclose{{}\left|\dot{~{}y}_{E}}\right|$ $\displaystyle\leq\mathopen{}\mathclose{{}\left(L_{\nabla E}\sup_{~{}x\in K} \mathopen{}\mathclose{{}\left|~{}L}\right|+L_{\tilde{~{}L}}\sup_{~{}x\in K} \mathopen{}\mathclose{{}\left|\nabla E}\right|}\right)\mathopen{}\mathclose{{} \left|~{}y}\right|+\varepsilon\mathopen{}\mathclose{{}\left(\sup_{~{}x\in K} \mathopen{}\mathclose{{}\left|\tilde{~{}L}}\right|+\sup_{~{}x\in K}\mathopen{} \mathclose{{}\left|\nabla E}\right|}\right)=:a_{E}\mathopen{}\mathclose{{} \left|~{}y}\right|+\varepsilon\,b_{E}.$ $$

Similarly, by adding and subtracting $~{}M(~{}x)\nabla S\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right),\tilde{~{
}M}(~{}x)\nabla S\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right),\tilde{~{}
M}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)\nabla S\mathopen{}
\mathclose{{}\left(\tilde{~{}x}}\right)$, it follows that

$$ $\displaystyle\dot{~{}y}_{S}$ $\displaystyle=~{}M(~{}x)\mathopen{}\mathclose{{}\left(\nabla S(~{}x)-\nabla S \mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)}\right)+\mathopen{} \mathclose{{}\left(~{}M(~{}x)-\tilde{~{}M}(~{}x)}\right)\nabla S\mathopen{} \mathclose{{}\left(\tilde{~{}x}}\right)$ $\displaystyle\qquad+\mathopen{}\mathclose{{}\left(\tilde{~{}M}(~{}x)-\tilde{~{ }M}\mathopen{}\mathclose{{}\left(\tilde{~{}x}}\right)}\right)\nabla S\mathopen {}\mathclose{{}\left(\tilde{~{}x}}\right)+\tilde{~{}M}\mathopen{}\mathclose{{} \left(\tilde{~{}x}}\right)\mathopen{}\mathclose{{}\left(\nabla S\mathopen{} \mathclose{{}\left(\tilde{~{}x}}\right)-\nabla\tilde{S}\mathopen{}\mathclose{{ }\left(\tilde{~{}x}}\right)}\right).$ $$

By Proposition 3.7, there exists a two-layer network $\tilde{~{}M}$ with one continuous derivative such that $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|~{}M-\tilde{~{}M}}\right|<\varepsilon$, with $\tilde{~{}M}$ Lipschitz continuous for the same reason as before. It follows from this and $\sup_{~{}x\in K}\mathopen{}\mathclose{{}\left|\nabla S-\nabla\tilde{S}}\right|<\varepsilon$ that

$$ $\displaystyle\mathopen{}\mathclose{{}\left|\dot{~{}y}_{S}}\right|$ $\displaystyle\leq\mathopen{}\mathclose{{}\left(L_{\nabla S}\sup_{~{}x\in K} \mathopen{}\mathclose{{}\left|~{}M}\right|+L_{\tilde{~{}M}}\sup_{~{}x\in K} \mathopen{}\mathclose{{}\left|\nabla S}\right|}\right)\mathopen{}\mathclose{{} \left|~{}y}\right|+\varepsilon\mathopen{}\mathclose{{}\left(\sup_{~{}x\in K} \mathopen{}\mathclose{{}\left|\tilde{~{}M}}\right|+\sup_{~{}x\in K}\mathopen{} \mathclose{{}\left|\nabla S}\right|}\right)=:a_{S}\mathopen{}\mathclose{{} \left|~{}y}\right|+\varepsilon\,b_{S}.$ $$

Now, recall that $\partial_{t}\mathopen{}\mathclose{{}\left|~{}y}\right|=\mathopen{}\mathclose{{
}\left|~{}y}\right|^{-1}\mathopen{}\mathclose{{}\left(\dot{~{}y}\cdot~{}y}
\right)\leq\mathopen{}\mathclose{{}\left|\dot{~{}y}}\right|$ by Cauchy-Schwarz, and therefore the time derivative of $\mathopen{}\mathclose{{}\left|~{}y}\right|$ is bounded by

$$ $\displaystyle\partial_{t}\mathopen{}\mathclose{{}\left|~{}y}\right|$ $\displaystyle\leq\mathopen{}\mathclose{{}\left|\dot{~{}y}_{E}}\right|+ \mathopen{}\mathclose{{}\left|\dot{~{}y}_{S}}\right|=\mathopen{}\mathclose{{} \left(a_{E}+a_{S}}\right)\mathopen{}\mathclose{{}\left|~{}y}\right|+ \varepsilon\mathopen{}\mathclose{{}\left(b_{E}+b_{S}}\right)=:a\mathopen{} \mathclose{{}\left|~{}y}\right|+b.$ $$

This implies that $\partial_{t}\mathopen{}\mathclose{{}\left|~{}y}\right|-a\mathopen{}\mathclose{
{}\left|~{}y}\right|\leq b$, so multiplying by the integrating factor $e^{-at}$ and integrating in time yields

$$ $\displaystyle\mathopen{}\mathclose{{}\left|~{}y(t)}\right|\leq\varepsilon b \int_{0}^{t}e^{a\mathopen{}\mathclose{{}\left(t-\tau}\right)}\,d\tau= \varepsilon\frac{b}{a}\mathopen{}\mathclose{{}\left(e^{at}-1}\right),$ $$

where we used that $~{}y(0)=~{}0$ since the initial condition of the trajectories is shared. Therefore, the $L^{2}$ error in time can be approximated by

$$ $\displaystyle\mathopen{}\mathclose{{}\left\|~{}y}\right\|^{2}$ $\displaystyle=\int_{0}^{T}\mathopen{}\mathclose{{}\left|~{}y}\right|^{2}\,dt \leq\varepsilon^{2}\frac{b^{2}}{a^{2}}\mathopen{}\mathclose{{}\left(e^{2aT}-2e ^{aT}+T+1}\right),$ $$

establishing the conclusion.
∎

###### Proof of Lemma 3.2 .

###### Proof of Theorem 3.4 .

###### Lemma A.1 .

###### Proof.

###### Lemma A.2 .

###### Proof.

###### Proof of Proposition 3.7 .

###### Proof of Theorem 3.9 .

## Appendix B Experimental and Implementation Details

This Appendix records additional details related to the numerical experiments in Section 5. For each benchmark problem, a set of trajectories is manufactured given initial conditions by simulating ODEs with known metriplectic structure. For the experiments in Table 1, only the observable variables are used to construct datasets, since entropic information is assumed to be unknown. Algorithm 2 summarizes the training of the dynamics models used for comparison with NMS.

Figure: Algorithm 2 Training dynamics models

For each compared method, integrating the ODEs is done via the Dormand–Prince method (dopri5) [] with relative tolerance $10^{-7}$ and absolute tolerance $10^{-9}$. The loss is evaluated by measuring the discrepancy between the ground truth observable states $~{}x^{\text{o}}$ and the approximate observable states $\tilde{~{}x}^{\text{o}}$ in the mean absolute error (MAE) metric. The model parameters $\Theta$ (i.e., the weights and biases) are updated by using Adamax [] with an initial learning rate of 0.01. The number of training steps is set as 30,000, and the model parameters resulting in the best performance for the validation set are chosen for testing.
Specific information related to the experiments in Section 5 is given in the subsections below.

For generating the results reported in Table 1, we implemented the proposed algorithm in Python 3.9.12 and PyTorch 2.0.0. Other required information is provided with the accompanying code. All experiments are conducted on Apple M2 Max chips with 96 GB memory. To provide the mean and the standard deviation, experiments are repeated three times with varying random seeds for all considered methods.

### B.1 Two gas containers

As mentioned in the body, the two gas container (TGC) problem tests models’ predictive capability (i.e., extrapolation in time). To this end, one simulated trajectory is obtained by solving an IVP with a known TGC system and an initial condition, and the trajectory of the observable variables is split into three subsequences, [0, $t_{\text{train}}$], ($t_{\text{train}},t_{\text{val}}]$, and ($t_{\text{val}},t_{\text{test}}]$ for training, validation, and test with $0<t_{\text{train}}<t_{\text{val}}<t_{\text{test}}$.

In the experiment, a sequence of 100,000 timesteps is generated using the Runge–Kutta 4th-order (RK4) time integrator with a step size 0.001. The initial condition is given as $~{}x=(1,2,103.2874,103.2874)$ following []. The training/validation/test split is defined by $t_{\text{train}}=20$, $t_{\text{val}}=30$, and $t_{\text{test}}=100$. For a fair comparison, all considered models are set to have a similar number of model parameters, $\sim$2,000.
The specifications of the network architectures are:

- •
NMS: The total number of model parameters is 1959. The functions $~{}A_{\mathrm{tri}},~{}B,~{}K_{\mathrm{chol}},E,S$ are parameterized as MLPs with the Tanh nonlinear activation function. The MLPs parameterizing $~{}A_{\mathrm{tri}},~{}B,~{}K_{\mathrm{chol}},E$ are specified as 1 hidden layer with 10 neurons, and the on parameterizing $S$ is specified as 3 hidden layers with 25 neurons.
- •
NODE: The total number of model parameters is 2179. The black-box NODE is parameterized as an MLP with the Tanh nonlinear activation function, 4 hidden layers and 25 neurons.
- •
SPNN: The total number of model parameters is 1954. The functions $E$ and $S$ are parameterized as MLPs with the Tanh nonlinear activation function; each MLP is specified as 3 hidden layers and 20 neurons. The two 2-tensors defining $~{}L$ and $~{}M$ are defined as learnable $3\times 3$ matrices.
- •
GNODE: The total number of model parameters is 2343. The functions $E$ and $S$ are parameterized as MLPs with the Tanh nonlinear activaton function; each MLP is specified as 2 hidden layers and 30 neurons. The matrices and 3-tensors required to learn $~{}L$ and $~{}M$ are defined as learnable $3\times 3$ matrices and $3\times 3\times 3$ tensor.
- •
GFINN: The total number of model parameters is 2065. The functions $E$ and $S$ are parameterized as MLPs with Tanh nonlinear activation function; each MLP is specified as 2 hidden layers and 20 neurons. The matrices to required to learn $~{}L$ and $~{}M$ are defined as $K$ learnable $3\times 3$ matrices, where $K$ is set to 2.

### B.2 Thermoelastic double pendulum

The equations of motion in this case are given for $1\leq i\leq 2$ as

$$ $\displaystyle\dot{~{}q}_{i}$ $\displaystyle=\frac{~{}p_{i}}{m_{i}},\quad\dot{~{}p}_{i}=-\partial_{~{}q_{i}} \mathopen{}\mathclose{{}\left(E_{1}(~{}x)+E_{2}(~{}x)}\right),\quad\dot{S}_{1} =\kappa\mathopen{}\mathclose{{}\left(T_{1}^{-1}T_{2}-1}\right),\quad\dot{S}_{2 }=\kappa\mathopen{}\mathclose{{}\left(T_{1}T_{2}^{-1}-1}\right),$ $$

where $\kappa>0$ is a thermal conductivity constant (set to 1), $m_{i}$ is the mass of the $i^{\mathrm{th}}$ spring (also set to 1) and $T_{i}=\partial_{S_{i}}E_{i}$ is its absolute temperature. In this case, $~{}q_{i},~{}p_{i}\in\mathbb{R}^{2}$ represent the position and momentum of the $i^{\mathrm{th}}$ mass, while $S_{i}$ represents the entropy of the $i^{\mathrm{th}}$ pendulum. As before, the total entropy $S(~{}x)=S_{1}+S_{2}$ is the sum of the entropies of the two springs, while defining the internal energies

$$ $\displaystyle E_{i}(~{}x)$ $\displaystyle=\frac{1}{2}\mathopen{}\mathclose{{}\left(\ln\lambda_{i}}\right)^ {2}+\ln\lambda_{i}+e^{S_{i}-\ln\lambda_{i}}-1,\quad\lambda_{1}=\mathopen{} \mathclose{{}\left|~{}q_{i}}\right|,\quad\lambda_{2}=\mathopen{}\mathclose{{} \left|~{}q_{2}-~{}q_{1}}\right|,$ $$

leads to the total energy $E(~{}x)=(1/2m_{1})\mathopen{}\mathclose{{}\left|~{}p_{1}}\right|^{2}+(1/2m_{2}
)\mathopen{}\mathclose{{}\left|~{}p_{2}}\right|^{2}+E_{1}(~{}x)+E_{2}(~{}x).$

The thermoelastic double pendulum experiment tests model prediction across initial conditions. In this case, 100 trajectories are generated by varying initial conditions that are randomly sampled from [0.1,1.1] $\times$ [-0.1,0.1] $\times$ [2.1, 2.3] $\times$ [-0.1,0.1] $\times$ [-1.9,2.1] $\times$ [0.9,1.1] $\times$ [-0.1, 0.1] $\times$ [0.9,1.1] $\times$ [0.1,0.3] $\subset\mathbb{R}^{10}$. Each trajectory is obtained from the numerical integration of the ODEs using an RK4 time integrator with step size 0.02 and the final time $T=40$, resulting in the trajectories of length 2,000. The resulting 100 trajectories are split into 80/10/10 for training/validation/test sets. For a fair comparison, all considered models are again set to have similar number of model parameters, $\sim$2,000. The specifications of the network architectures are:

- •
NMS: The total number of model parameters is 2201. The functions $~{}A,~{}B,~{}K,E,S$ are parameterized as MLPs with the Tanh nonlinear activation function. The MLPs parameterizing are specified as 1 hidden layer with 15 neurons.
- •
NODE: The total number of model parameters is 2005. The black-box NODE is parameterized as an MLP with the Tanh nonlinear activation function, 2 hidden layers and 35 neurons.
- •
SPNN: The total number of model parameters is 2362. The functions $E$ and $S$ are parameterized as MLPs with the Tanh nonlinear activation function; each MLP is specified as 3 hidden layers and 20 neurons. The two 2-tensors defining $~{}L$ and $~{}M$ are defined as learnable $3\times 3$ matrices.
- •
GNODE: The total number of model parameters is 2151. The functions $E$ and $S$ are parameterized as MLPs with the Tanh nonlinear activaton function; each MLP is specified as 2 hidden layers and 15 neurons. The matrices and 3-tensors required to learn $~{}L$ and $~{}M$ are defined as learnable $3\times 3$ matrices and $3\times 3\times 3$ tensor.
- •
GFINN: The total number of model parameters is 2180. The functions $E$ and $S$ are parameterized as MLPs with Tanh nonlinear activation function; each MLP is specified as 2 hidden layers and 15 neurons. The matrices to required to learn $~{}L$ and $~{}M$ are defined as $K$ learnable $3\times 3$ matrices, where $K$ is set to 2.

## Appendix C Additional experiment: Damped nonlinear oscillator

Consider a damped nonlinear oscillator of variable dimension with state $~{}x=\begin{pmatrix}~{}q&~{}p&S\end{pmatrix}^{\intercal}$, whose motion is governed by the metriplectic system

$$ $\displaystyle\dot{~{}q}=\frac{~{}p}{m},\quad\dot{~{}p}=k\sin{~{}q}-\gamma~{}p, \quad\dot{S}=\frac{\gamma\mathopen{}\mathclose{{}\left|~{}q}\right|^{2}}{mT}.$ $$

Here $~{}q,~{}p\in\mathbb{R}^{n}$ denote the position and momentum of the oscillator, $S$ is the entropy of a surrounding thermal bath, and the constant parameters $m,\gamma,T$ are the mass, damping rate, and (constant) temperature. This leads to the total energy $E(~{}x)=(1/2m)\mathopen{}\mathclose{{}\left|~{}p}\right|^{2}-k\cos{~{}q}+TS,$
which is readily seen to be constant along solutions $~{}x(t)$.

It is now verified that NMS can accurately and stably predict the dynamics of a nonlinear oscillator $~{}x=\begin{pmatrix}~{}q&~{}p&S\end{pmatrix}^{\intercal}$ in the case that $n=1,2$, both when the entropy $S$ is observable as well as when it is not. As before, the task considered is prediction in time, although all compared methods NODE, GNODE, and $\text{NMS}_{\text{known}}$ are now trained on full state information from the training interval, and test errors are computed over the full state $~{}x$ on the extrapolation interval $(t_{\mathrm{valid}},t_{\mathrm{test}}]$, which is $150\%$ longer than the training interval. In addition, another NMS model, $\text{NMS}_{\mathrm{diff}}$, was trained using only the partial state information $~{}x^{o}=\begin{pmatrix}~{}q,~{}p\end{pmatrix}^{\intercal}$and tested under the same conditions, with the initial guess for $~{}x^{u}$ generated as in Appendix E.
As can be seen in Table 2, NMS is more accurate than GNODE or NODE in both the 1-D and 2-D nonlinear oscillator experiments, improving on previous results by up to two orders of magnitude. Remarkably, NMS produces more accurate entropic dynamics even in the case where the entropic variable $S$ is unobserved during NMS training and observed during the training of other methods. This illustrates another advantage of the NMS approach: because of the reasonable initial data for $S$ produced by the diffusion model, the learned metriplectic system produced by NMS remains performant even when metriplectic governing equations are unknown and only partial state information is observed.

**Table 2: Experimental results for the benchmark problems with respect to MSE and MAE. The best scores are in boldface.**
|  | 1-D D.N.O. | T.G.C. | 2-D D.N.O. |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| MSE | MAE | MSE | MAE | MSE | MAE |  |
| $\textbf{NMS}_{\mathrm{diff}}$ | .0170 | .1132 | .0045 | .0548 | .0275 | .1456 |
| $\textbf{NMS}_{\mathrm{known}}$ | .0239 | .1011 | .0012 | .0276 | .0018 | .0357 |
| NODE | .0631 | .2236 | .0860 | .2551 | .0661 | .2096 |
| GNODE | .0607 | .1976 | .0071 | .0732 | .2272 | .4267 |

To describe the experimental setup precisely, data is collected from a single trajectory with initial condition as $~{}x=(~{}2,~{}0,0)$ following  []. The path is calculated at 180,000 steps with a time interval of 0.001, and is then split into training/validation/test sets as before using $t_{\text{train}}=60$, $t_{\text{val}}=90$ and $t_{\text{test}}=180$. Specifications of the networks used for the experiments in Table 2 are:

- •
NMS: The total number of parameters is 154. The number of layers for $~{}A_{\mathrm{tri}},~{}B,~{}K_{\mathrm{chol}},E,S$ is selected from {1,2,3} and the number of neurons per layer from {5,10,15}. The best hyperparameters are 1 hidden layer with 5 neurons for each network function.
- •
GNODE: The total number of model parameters is 203. The number of layers and number of neurons for each network is chosen from the same ranges as for NMS. The best hyperparameters are 1 layer with 10 neurons for each network function.
- •
NODE: The total number of model paramters is 3003. The NODE architecture is formed by stacking MLPs with Tanh activation functions. The number of blocks is chosen from {3,4,5} and the number of neurons of each MLP from {30,40,50}. The best hyperparameters are 4 and 30 for the number of blocks and number of neurons, respectively.

## Appendix D Scaling study

To compare the scalability of the proposed NMS architecture design with existing architectures, different realizations of GNODE, GFINN, and NMS are generated by varying the dimension of the state variables, $n=\{1,5,10,15,20,30,50\}$. The specifications of these models (i.e., hyperparameters) are set so that the number of model parameters is kept similar between each method for smaller values of $n$. For example, for $n=1,5$ the number of model parameters is $\sim$20,000 for each architecture. The results in Figure 2(a) confirm that GNODE scales cubically in $n$ while both GFINN and NMS scale quadratically. Note that only a constant scaling advantage of NMS over GFINN can be seen from this plot, since $r$ is fixed during this study.

It is also worthwhile to investigate the computational timings of these three models. Considering the same realizations of the models listed above, i.e., the model instances for varying $n=\{1,5,10,15,20,30,50\}$, 1,000 random samples of states $\{~{}x^{(i)}\}_{i=1}^{1,000}$ are generated. These samples are then fed to the dynamics function $~{}L(~{}x^{(i)})\nabla E(~{}x^{(i)})+~{}M(~{}x^{(i)})\nabla S(~{}x^{(i)})$ for $i=1,\ldots,1000$, and the computational wall time of the function evaluation via PyTorch’s profiler API is measured. The results of this procedure are displayed in Figure 2(b). Again, it is seen that the proposed NMSs require less computational resources than GNODEs and GFINNs.

Figure: (a) State dimension $n$ versus number of model parameters
Refer to caption: extracted/5627690/figs/n_vs_num_param_no_grid.png

Figure: (b) Model parameters versus wall time in microseconds
Refer to caption: extracted/5627690/figs/num_param_vs_wall_time_no_grid_loglog.png

## Appendix E Diffusion model for unobserved variables

Recent work in  [] suggests the benefits of performing time-series generation using a diffusion model. This Appendix describes how this technology is used to generate initial conditions for the unobserved NMS variables in the experiments corresponding to Table 2. More precisely, we describe how to train a conditional diffusion model which generates values for unobserved variables $~{}x^{u}$ given values for the observed variables $~{}x^{o}$.

#### Training and sampling:

Recall that diffusion models add noise with the following stochastic differential equation (SDE):

$$ $d\textbf{x}(t)=\textbf{f}(t,\textbf{x}(t))dt+g(t)d\textbf{w},\quad t\in[0,1],$ $$

where $\textbf{w}\in\mathbb{R}^{\dim(\textbf{x})}$ is a multi-dimensional Brownian motion, $\textbf{f}(t,\cdot):\mathbb{R}^{\dim(\textbf{x})}\rightarrow\mathbb{R}^{\dim(
\textbf{x})}$ is a vector-valued drift term, and $g:[0,1]\rightarrow\mathbb{R}$ is a scalar-valued diffusion function.

For the forward SDE, there exists a corresponding reverse SDE:

$$ $d\textbf{x}(t)=[\textbf{f}(t,\textbf{x}(t))-g^{2}(t)\nabla_{\textbf{x}(t)}{ \log p(\textbf{x}(t))}]dt+g(t)d\bar{\textbf{w}},$ $$

which produces samples from the initial distribution at $t=0$. This formula suggests that if the score function, $\nabla_{\textbf{x}(t)}{\log p(\textbf{x}(t))}$, is known, then real samples from the prior distribution $p(\textbf{x})\sim\mathcal{N}(\mu,\sigma^{2})$ can be recovered, where $\mu,\sigma$ vary depending on the forward SDE type.

In order for a model $M_{\theta}$ to learn the score function, it has to optimize the following loss:

$$ $L(\theta)=\mathbb{E}_{t}\{\lambda(t)\mathbb{E}_{\textbf{x}(t)}[\mathopen{} \mathclose{{}\left\|M_{\theta}(t,\textbf{x}(t))-{\nabla}_{\textbf{x}(t)}\log p (\textbf{x}(t))}\right\|_{2}^{2}]\},$ $$

where $t$ is uniformly sampled over $[0,1]$ with an appropriate weight function $\lambda(t):[0,1]\rightarrow\mathbb{R}$. However, using the above formula is computationally prohibitive.
Thanks to [], this loss can be substituted with the following denoising score matching loss:

$$ $\displaystyle L^{*}(\theta)$ $\displaystyle=\mathbb{E}_{t}\{\lambda(t)\mathbb{E}_{\textbf{x}(0)}\mathbb{E}_{ \textbf{x}(t)|\textbf{x}(0)}[\mathopen{}\mathclose{{}\left\|M_{\theta}(t, \textbf{x}(t))-{\nabla}_{\textbf{x}(t)}\log p(\textbf{x}(t)|\textbf{x}(0))} \right\|_{2}^{2}]\}.$ $$

Since score-based generative models use an affine drift term, the transition kernel $p(\textbf{x}(t)|\textbf{x}(0))$ follows a certain Gaussian distribution [], and therefore the gradient term ${\nabla}_{\textbf{x}(t)}\log p(\textbf{x}(t)|\textbf{x}(0))$ can be analytically calculated.

#### Experimental details

On the other hand, the present goal is to generate unobserved variables $~{}x^{u}$ given values for the observed variables $~{}x^{o}=(~{}q,~{}p)$, i.e., conditional generation. Therefore, our model has to learn the conditional score function, ${\nabla}_{~{}x^{u}(t)}\log p(~{}x^{u}(t)|~{}x^{o})$.
For example, in the damped nonlinear oscillator case, $S(t)$ is initialized as a perturbed $t\in[0,1]$, from which the model takes the concatenation of $~{}q,~{}p,S(t)$ as inputs and learns conditional the score function $\nabla_{S(t)}\log(S(t)|~{}q,~{}p)$.

For the experiments in Table 2, diffusion models are trained to generate $~{}x^{u}$ variables on three benchmark problems: the damped nonlinear oscillator, two gas containers, and thermolastic double pendulum. On each problem, representative parameters such as mass or thermal conductivity are varied, with the total number of cases denoted by $N$. Full trajectory data of length $T$ is then generated using a standard numerical integrator (e.g., dopri5), before it is evenly cut into $\lfloor{T/L}\rfloor$ pieces of length $L$. Let $V,U$ denote the total number of variables and the number of unobserved variables, respectively. It follows that the goal is to generate $U$ unobserved variables given $V-U$ observed ones, i.e., the objective is to generate data of shape $(NT/L,L,U)$ conditioned on data of shape $(NT/L,L,V-U)$. After the diffusion model has been trained for this task, the output data is reshaped into size $(N,T,U)$, which is used to initialize the NMS model. Note that the NODE and GNODE methods compared to NMS in Table 2 use full state information for their training, i.e., $~{}x^{u}=\varnothing$ in these cases, making it comparatively easier for these methods to learn system dynamics.

As in other diffusion models e.g. [], a U-net architecture is used, modifying 2-D convolutions to 1-D ones and following the detailed hyperparameters described in []. Note the following probability flow ODE seen in []:

$$ $\displaystyle d\textbf{x}(t)=\mathopen{}\mathclose{{}\left[\textbf{f}(t, \textbf{x}(t))-\frac{1}{2}g^{2}(t)\nabla_{\textbf{x}(t)}{\log p(\textbf{x}(t)) }}\right]dt,$ $$

Although models trained to mimic the probability flow ODE do not match the perofrmance of the forward SDE’s result in the image domain, the authors of [] observe that the probability flow ODE
outperforms the forward SDE in the time-series domain. Therefore, the probability flow ODE is used with the default hyperparameters of [].

## Peer Reviews

confidence: 2
contribution: 2
initial_rating: 5
presentation: 2
questions: See questions in the "Weaknesses" section above regarding physics background, training strategy, and alternative methods.
soundness: 3
strengths: The proposed neural network model for parameterizing the metriplectic system is novel to my knowledge, and it is well-motivated by the theoretical results on metriplectic operators and approximation errors presented in Section 3.2 and 3.4. Compared to several methods from prior literature, the proposed NMS method holds an advantage in terms of model size as well as empirical performances. The writing of the paper is clear overall.
summary: In this work, the authors present a neural network model (NMS) for learning the dynamics of metriplectic systems from trajectory data. It is based on a parameterization of certain scalar- and matrix-valued fields using neural networks that differs from prior works that enforce hard constraints on the degeneracy condition of metriplectic systems and results in a model size that scales quadratically in the system dimenion. The proposed method is demonstrated empirically to perform well in the learning of two low-dimensional metriplectic systems compared to prior methods.
weaknesses: **Physics background and relation to Hamiltonian systems:** Having relatively little prior knowledge about metriplectic systems, I would appreciate an expanded introduction of its physical motivation and some concrete examples for illustrating the general governing equation on the top of Page 2. For example, what do L and M look like in the two examples used for the numerical systems, or in general Hamiltonian systems -- will L be the "J" matrix in Hamiltonian systems (padded with zeros for the extra dimensions) and lose its independence on the state x? 

**On the decoupled block-wise structure:** The authors mentioned prior works such as Ruiz et al. (2021) and Xu et al. (2022 & 2023) which proposed to parameterize metriplectic systems assuming a decoupled block-wise structure and discussed their inability to express general metriplectic dynamics. But I would appreciate some examples of these more general metriplectic systems encountered in practice. In particular, do the two systems studied in Section 5 admit the decoupled block-wise structure? If so, it would be reasonable to expect that those more restrictive methods are also tested in the experiments as baselines.

**Comparison with Hamiltonian learning:** If one ignores the entropy states and focuses only on the observable states (positions and momenta), it looks like the two examples in Section 5 can just be learned as Hamiltonian systems. In that case, does the NMS method reduce to e.g. the Hamiltonian Neural Network (HNN) from [1]? If not, perhaps HNN should also be added as a baseline method to compare against in Table 2.

**Choices of the initial unobserved states:** Regarding the initial unobserved states in batch-wise training. The authors mentioned two interesting strategies to handle the missing initial unobserved entropy states (first question: which one was used in the experiments whose results are reported in the main text?), one of which is to assume that the entropy values increase linearly in time. Is there a justification behind this? (A further question seems to be: are the entropy states and their dynamics uniquely determined by the observable states?) Besides the two strategies considered by the authors, the initial state optimization proposed in [2] for learning Hamiltonian systems might also be an alternative to consider.

**Test systems are low-dimensional:** Another limitation, as acknowledged by the authors, is that the empirical performance of NMS on larger-scale, realistic metriplectic systems has not yet been demonstrated.

A minor issue: misplaced parentheses for citations on Page 1 (probably due to mixing up \citep with \citet).

References:

[1] Greydanus et al., "Hamiltonian Neural Networks", NeurIPS 2019.

[2] Chen et al., "Symplectic Recurrent Neural Networks.", ICLR 2020.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
questions: Besides lower complexity the method achieves better loss values compared with all the other methods. Why is that even on the simple two-gas problem the other methods, say GNODE, do not achieve the same accuracy even with full state information? Is it due to a restrictive parameterization assumption?  Some intuition would be useful.
soundness: 3
strengths: The proposed approach appears to be novel and more general than prior work and allows for all metriplectic data to be approximated simultaneously. Prior work assumes special forms to satisfy the constraints on entropy and energy. This allows the method to model a greater class of systems. The modeling assumptions are also quite mild and only require non-zero gradients for energy and entropy.

The claims and objectives of the paper are clear. The literature review is quite extensive and the experiments treat a number of prior baseline and validate the claims of the paper of better generalization at lower complexity.
summary: Metricplectic systems model systems that satisfy energy conservation as well as entropy constraints and can used to model thermodynamic and other system that require such constraints. In the context of machine learning it is of interest to learn such systems from data.
Prior methods like GNODE exchange the problem of enforcing degeneracy constraints with the problem of enforcing symmetry constraints which underdetermines the problem and at the same time have a redundant parameterization of the problem which leads to high (cubic) complexity.
The proposed method exploits structure in the tensor fields to reduce the number of parameters.
The paper demonstrates structure in the degeneracy constraints beyond what can be captured by symmetry constraints leading to a lower parameter parameterization.
A further result shows that the proposed formulation universally approximates metricplectic systems that are non-degenerate and shows a generalization result.
weaknesses: The paper is difficult to read and there is not a lot of intuition to understand the source of the complexity reduction in Lemma 3.2 and Theorem 3.4 for a non-expert.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 2
questions: 0. How original is the approach to using exterior algebra to parameterize hard-constrained structures? The reviewer believes some papers that deal with structure preservation using exterior algebra exist in the literature. If the authors also know such papers, the authors could have remarked on that.

1. In Remark 3.5, it is claimed that "the proposed parameterizations for L, M are not one-to-one but properly contain the set of valid nondegenerate metriplectic systems". Are the authors trying to say that the proposed parameterizations are not injective but subjective onto the set of possible nondegenerate metriplectic systems? I think so from the proof. The author could clarify what they mean by "properly contain".

2. Again in Remark 3.5, the authors said that the Jacobi identity is not enforced in the algorithm, which causes the parameterization to be not one-to-one. It is not clear to the reviewer how are these related. From the reviewer's viewpoint, the parameterization is not one-to-one because the construction is via orthogonal projection, then of course there will be more than one parameterizations that give the same metriplectic system. However, it is unclear how this is linked to the Jacobi identity not being enforced. The author should clarify this point. Besides, in the last sentence of Remark 3.5, it is said that the structure and energy conservation cannot be simultaneously preserved, which, as far as the reviewer knows, applies in the case of symplectic integrators, but here the context is not symplectic integrators. The authors could clarify this point.

3. In page 6, it is claimed that "the exterior algebraic expressions in Lemma 3.2 require less redundant operations than the corresponding metricized expressions from Theorem 3.4, and therefore the expressions from Lemma 3.2 are used when implementing NMS". The author should clarify how the construction in Lemma 3.2 requires less redundant operations.

4. On top of page 8 (lines 381-382), the first strategy to deal with unobserved states is to assume a line between the all 0 vector to the all 1 vector. Can this be justified?

5. In Algorithm 1, the input, xs = x(ts, μs), but what is μs? It seems that it is not introduced.

6. There is a typo "Lemmata" in line 641.

7. The algorithm in the end still needs structure-preserving integration, which means for each system under consideration, a structure-preserving integrator needs to be applied. In the numerical experiment, how are the integrators chosen?
soundness: 3
strengths: 1. The subject being investigated is of significant importance in AI for Science. The topic is closely related to the so-called structure-preserving machine learning. The paper is in general well written.

2. The idea of the proposed parameterization is novel. The construction of the parameterizations that satisfy the hard constraints is obtained through exterior algebra computations. The construction can be seen as obtained by orthogonal projection, which is quite natural. In Theorem 3.4, the statement is "if and only if", which is a nice result for machine learning purposes. There are universality analysis and growth rate of error analysis, which are theoretically convincing.

3. The numerical experiments are carried out quite extensively, together with comparisons with other machine learning schemes in the literature.
summary: The paper presents a new parameterization scheme for learning metriplectic systems. The parameterization is constructed based on computations using properties of the exterior algebra. The constraints can be seen as hard constraints instead of soft constraints. Proof of the construction and proof of universal approximation are provided as well as proof of growth rate of error with respect to time. An algorithm has been proposed, and numerical experiments are carried out to justify the effectiveness of the proposed scheme.
weaknesses: The main weakness of the paper is that, some of the remarks and claims are not clearly explained. The reviewer will state the details in the "Question" session.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: 1. While it is evident that the proposed parameterization requires fewer learnable scalar functions, this alone does not directly imply superior performance. Could the paper provide an explanation of why the proposed approach achieves higher accuracy?

2. The paper claims that the core advantage of the proposed method is its efficiency. Could the paper provide evidence demonstrating that the method requires less computational time or other resources, such as memory?

3. Could the paper provide that specific formula of metric use to compare the performance?

4. As the true governing functions (the right hand side of the ODE) are known, could the paper show the results of learning these governing functions?
soundness: 4
strengths: 1.  The emphasis on energy conservation and entropy stability is crucial for modeling realistic physical systems. The authors demonstrate that their approach maintains these properties, which is essential for ensuring that the learned models are physically meaningful and applicable to real-world scenarios

2. The proposed parameterization of metriplectic systems requires less learnable scalar functions than existing methods. 

3. The paper is well-structured and clearly written, making complex concepts accessible to a broader audience. The examples and comparisons with previous methods help to illustrate the advantages of the proposed approach effectively.
summary: This paper introduces a novel framework for learning metriplectic dynamics, which describe systems that conserve energy while generating entropy. The key contribution of the proposed approach lies in its design, which ensures both energy conservation and entropy stability. Mathematically, the paper presents a parameterization method that employs neural networks to model four functions, $L$, $E$, $M$, and $S$, such that $L \nabla S = M \nabla E = 0$. Here, $L$ is an antisymmetric matrix-valued function, $M$ is a symmetric matrix-valued function, and $E$ and $S$ are scalar functions.

Under the assumption that $\nabla E,\nabla S\neq 0$, the proposed parameterization of metriplectic systems requires less learnable scalar functions than existing methods. Numerical results show that  the proposed approach exhibits superior accuracy that existing methods for  learning metriplectic dynamics from data.
weaknesses: 1. Although the paper presents empirical results demonstrating the model's performance, the range of examples may be limited. Specifically, in both examples, the entropy is given by $S = s_1 + s_2$, resulting in a constant gradient $\nabla S = 0 $. This contradicts the core assumption in the paper that $ \nabla S \neq 0 $. 

2. The requirement that the metriplectic system being approximated is nondegenerate—i.e., the gradients of energy and entropy must not vanish ($\nabla E, \nabla S \neq 0$)—may limit the applicability of the method. Although this is claimed to be a mild condition, both of the examples considered in the paper contradict this core assumption. In addition to the entropy $S$ having a zero gradient ($\nabla S = 0$), the energy in the investigated examples is also degenerate. From the reviewers' understanding, steady-state is a fundamental concept in physical systems, characterized by the energy reaching an extremum, meaning the energy gradient is zero ($\nabla E = 0$). The reviewers are not aware of any physical systems with non-zero energy gradients. To better demonstrate the applicability of the proposed method, the reviewers recommend that the paper provide a detailed discussion of systems that satisfy the key non-degeneracy assumption. Additionally, it would be beneficial to include examples that adhere to this assumption to validate the effectiveness of the algorithm.

Meta-review: This paper proposes a method for learning the dynamics of metriplectic systems from trajectory data. It hardwires structural information in the learning model using geometric techniques, so that the physical constraint is strictly satisfied. Reviewers and I found the method to be novel. Although reviewers (e.g., gp5d and QHfC) initially expressed some concerns about the motivation, presentation and applicability, many of these concerns were alleviated during the discussions. Overall, I feel the pros outweigh the cons and therefore recommend acceptance.

Decision: Accept (Poster)