# Learning of discrete models of variational PDEs from data

Christian Offen ${\href https://orcid.org/0000-0002-5940-8057}$
<EMAIL_ADDRESS>https://www.uni-paderborn.de/en/person/85279
Sina Ober-Blöbaum  Paderborn University, Department of Mathematics, Warburger
Str. 100, 33098 Paderborn, Germany

###### Abstract

We show how to learn discrete field theories from observational data of fields
on a space-time lattice. For this, we train a neural network model of a
discrete Lagrangian density such that the discrete Euler–Lagrange equations
are consistent with the given training data. We, thus, obtain a structure-
preserving machine learning architecture. Lagrangian densities are not
uniquely defined by the solutions of a field theory. We introduce a technique
to derive regularisers for the training process which optimise numerical
regularity of the discrete field theory. Minimisation of the regularisers
guarantees that close to the training data the discrete field theory behaves
robust and efficient when used in numerical simulations. Further, we show how
to identify structurally simple solutions of the underlying continuous field
theory such as travelling waves. This is possible even when travelling waves
are not present in the training data. This is compared to data-driven model
order reduction based approaches, which struggle to identify suitable latent
spaces containing structurally simple solutions when these are not present in
the training data. Ideas are demonstrated on examples based on the wave
equation and the Schrödinger equation.

††preprint: AIP/123-QED

> To learn practical models of dynamical systems from data, prior geometric
> knowledge and aspects of numerical integration theory can be exploited in
> the design of the data-driven architectures. The article shows how to
> utilise numerical analysis informed neural networks and discrete field
> theories to model dynamical systems governed by first principles.

##  I Introduction

The identification of models of dynamical systems from observational data is
an important task in machine learning: the identified models can be used to
interpolate data points, predict motions, develop control strategies, or
analyse long-term behaviourNarendra and Parthasarathy (1990); Han  _et al._
(2019). Additionally, physical properties of the system can be discovered such
as symmetries, conservation lawsDierkes  _et al._ (2023); Lishkova  _et al._
(2023), stress tensorsWu, Xiao, and Paterson (2018), for instance. Moreover,
data-driven models can be used as surrogates for analytic models when
computations with explicit analytical models are prohibitively expensive, such
as in molecular dynamicsKadupitiya  _et al._ (2020).

To improve qualitative aspects of data-driven models, prior knowledge about
physical aspects of the dynamics, such as symmetries, conservation laws,
Hamiltonian or variational structure can be taken into account. The approach
is referred to as physics informed machine learning. (Notice, however, that
the term physics informed neural networks has also been used to describe
neural network approximations of solutions to partial differential equations
that represent physical laws. See Karniadakis et alKarniadakis  _et al._
(2021) for a review article.)

Variational principles are at the hart of physical theories and constitute
first principles from which dynamical models can be derived. In such dynamical
systems motions or solutions to the field equations constitute stationary
points of an action functional. For instance, classical mechanical systems,
electromagnetic field theories, fluid dynamics, as well as quantum dynamical
systems are governed by variational principles. Presence of variational
structure in a dynamical system is related to many profound laws in physics
such as the correspondence of symmetries and conserved quantities (Noether’s
theorem) Marsden and Ratiu (1999); Mansfield (2010).

To embed variational structure into machine learned models of dynamical
systems, Greydanus et al propose to learn a model of a variational principle
from data, called Lagrangian neural networks (LNNs) Cranmer  _et al._ (2020).
This needs to be contrasted to learning of a model of a flow map or a vector
field of a dynamical system without enforcing physical properties of the
system Hamzi and Owhadi (2021). Similar ideas and extensions to LNN have been
developed, such as combinations with model order reduction techniques Allen-
Blanchette  _et al._ (2020) or the inclusion of external forces Lutter,
Ritter, and Peters (2019).

Other approaches focus on exploiting Hamiltonian structure for learning models
of dynamical systems. While Hamiltonian and variational structure are
equivalent from a theoretical perspective (under non-degeneracy conditions),
data-driven architectures for learning Hamiltonian systems, such as
Hamiltonian Neural NetworksGreydanus, Dzamba, and Yosinski (2019) or
Symplectic Neural Networks Jin  _et al._ (2020), differ from those
architectures that learn models of action functionals as they require
additional observation of momentum data rather than just position data (and
possibly derivatives), or at least some prior knowledge of the symplectic
phase space structure. However, much of the symplectic structure can be
learned from data as well Bertalan  _et al._ (2019); Chen  _et al._ (2023).

Instead of learning a model of a continuous action functional to describe a
dynamical system with variational structure, QinQin (2020) proposes to learn a
discrete variational principle instead. Discrete field theories can be trained
on discrete observations of fields over a space-time lattice, while the
training of continuous theories such as LNNCranmer  _et al._ (2020) requires
observations of velocities and acceleration data. Moreover, solutions of
discrete field theories can be computed without further discretisation by
numerical integrators such that discretisation errors are avoided.

Figure 1: Left: Motions of a data-driven model (solid lines) for the
mathematical pendulum match motions of an analytic model (dotted). Centre:
Numerically computed motions of the data-driven model are very inaccurate.
Right: Applying the same numerical integrator to the analytic model yield much
more accurate results. (Figure taken from (Ober-Blöbaum and Offen, 2023, Fig.
6). See Ober-Blöbaum and Offen (2023) for details.)

However, the authors have demonstrated that numerical discretisation errors of
learned continuous variational models of dynamical systems can be surprisingly
bigOber-Blöbaum and Offen (2023): the plot to the left of Fig. 1 shows that
the exact trajectories of the learned dynamical system of the mathematical
pendulum coincide with the exact trajectories of an analytic reference.
However, an application of a standard variational numerical integrator
(variational midpoint ruleMarsden and West (2001)) to the learned model
performs poorly (centre) while an application of the same integrator with the
same step-size shows only small numerical errors when applied to the analytic
reference (right). The reason is that variational principles are not uniquely
determined by a system’s motions. In machine learning, care is required to
identify a model of a variational principle that is not only consistent with
the motion data, but also suitable for numerical computations. This will be a
central aspect of this article.

While learning a discrete field theory avoids discretisation errors in
numerical computation, additional regularisation is essential to avoid
learning degenerate principles for which discrete field equations become
unsolvable or badly conditioned. While the authors have discussed a geometric
regularisation strategy applying to models based on Gaussian Processes or
kernel methodsOber-Blöbaum and Offen (2023), different regularisation
strategies are required when artificial neural networks are used. For
variational dynamical systems, for which the underlying continuous system is
an ordinary differential equation, we show first approaches in Lishkova et
alLishkova  _et al._ (2023). The present article extends the authors’
conference contributionOffen and Ober-Blöbaum (2023), in which machine
learning architectures for discrete field theories are developed whose
underlying dynamical system is a partial differential equation.

Learning of discrete field theories governing infinite-dimensional dynamical
systems directly from data needs to be contrasted with approaches based on
model order reduction (MOR). MOR-based techniques learn a data-driven
reduction map of a spatially discretised system. Then, a full order model
(FOM) is projected to the latent space to obtain a reduced order model (ROM).
If no FOM is available, the learned reduction map is used to project dynamical
data to the latent space. Subsequently, the ROM is learned, i.e. a data-driven
model is fitted to describe the temporal dynamics of the system on the latent
space. Several approaches have been developed to preserve geometric structure
in this contextAllen-Blanchette  _et al._ (2020); Mason  _et al._ (2022);
Buchfink, Glas, and Haasdonk (2023); Carlberg, Tuminaro, and Boggs (2015);
Sharma  _et al._ (2023); Sharma and Kramer (2022); Sharma, Wang, and Kramer
(2022); Tyranowski and Kraus (2022).

In contrast to the MOR-based techniques, we learn the discrete field theory
directly without a reduction step. Assuming that the true dynamical system can
be described by variational partial differential equations, the true,
continuous field theory is local in all variables (space and time) in the
following sense: whether a function $u$ constitutes a solution to the field
theory can be tested locally, by checking whether $u$ fulfils a partial
differential equation at any point in the space-time domain. MOR-based data-
driven models do not preserve the locality of the field theory in all
variables. When reduction is applied to the spatial directions, the resulting
ROM is local in time but not in space. Therefore, unphysical global effects in
the spatial directions are not excluded by design. In contrast, discrete field
theories are local in all variables (space and time) by construction. We
exploit this locality to obtain a reduction of computational complexity, which
would otherwise require MOR.

Moreover, we demonstrate that avoiding a projection step can be advantageous
if, for instance, highly symmetric solutions of the model, such as travelling
waves, need to be identified in the data-driven model. We develop a method to
detect and preserve solutions of the discrete field theory that correspond to
highly symmetric solutions in the underlying continuous dynamical system. As
the mesh of a discrete field theory is typically not compatible with the
symmetries of the continuous system, we use a suitable symmetric ansatz of the
discrete action functional based on Palais’ principle of symmetric criticality
Palais (1979); Torre (2011) to design machine learning architectures for
discrete field theories with travelling waves.

The main novelties of the article are

  * •

the systematic design of numerical analysis informed regularisation strategies
for machine learning models of discrete field theories,

  * •

the preservation and detection of highly symmetric solutions in data-driven
discrete field theories using data-driven architectures based on Palais’
principle of symmetric criticality and comparison to data-driven model order
reduction based techniques,

  * •

extension of the approaches first introduced in our conference
contributionOffen and Ober-Blöbaum (2023) to different stencils and to
degenerate Lagrangians that are linear in velocities.

The remainder of the article is structured as follows: Section II reviews the
concept of variational principles in the continuous and discrete setting. The
theory is illustrated on our main examples, the wave equation and the
Schrödinger equation. Section III provides theoretical considerations
motivating the regularisers that we will employ when training our neural
network models. Section IV introduces the machine learning set-up to learn a
discrete field theory from data. Section IV provides theoretical background on
the notion of travelling waves and highly symmetric solutions in discrete
field theories and describes their variational structure. Finally, Section VI
provides numerical experiments based on the wave and Schrödinger equation. The
article concludes with a summery and future work section (Section VII).

##  II Continuous and discrete variational principles and field theories

We recall briefly the concept of variational principles and their
discretisations. We illustrate the theory on the wave equation and the
Schrödinger equation, on which we will test our data-driven framework. For
more detailed expositions, we refer to the literature for an introduction to
variational calculusOlver (1986); Basdevant (2007); Mansfield (2010) and
discrete mechanicsMarsden and West (2001).

###  II.1 Continuous variational principles

Variational principles play a fundamental role in the derivation and analysis
of field theories in physicsBasdevant (2007): motions $u\colon
X\to\mathbb{R}^{d}$ are described as stationary points of an action functional

$S(u)=\int_{X}L\left(x,u(x),u_{x_{0}}(x),\ldots,u_{x_{n}}(x)\right)\mathrm{d}x.$
(1)

Here the domain of definition of $u$, $X$, refers to a suitable manifold (with
boundary) such as
$X=[x^{0}_{0},x^{N}_{0}]\times\ldots\times[x^{0}_{n},x^{N}_{n}]$ with
coordinates $x=(x_{0},x_{1},\ldots,x_{n})$. Often, $X$ represents a space-time
domain. The expression $u_{x_{k}}$ denotes the partial derivative
$\frac{\partial u}{\partial x_{k}}$. More precisely, $u$ constitutes a
solution of the field theory, if the action functional $S$ is stationary at
$u$ with respect to variations $\delta u\colon X\to\mathbb{R}^{d}$ vanishing
at the boundary of $X$. The function $L$ is called Lagrangian density. As it
depends on the 1-jet of $u$ only and not on higher derivatives, a field theory
described by an action functional of the from (1) is called 1st order field
theory. Stationary points of the action functional are solutions to the Euler-
Lagrange equation

$0=\mathrm{EL}(L)=\frac{\partial L}{\partial
u}-\sum_{j=0}^{n}\frac{\mathrm{d}}{\mathrm{d}x_{j}}\frac{\partial L}{\partial
u_{x_{j}}},$ (2)

which constitutes a partial differential equation of 2nd order, unless $L$ is
degenerate.

######  Remark 1 (Gauge freedom).

For a continuously differentiable function $F=(F^{1},\ldots,F^{n})$ in
$(x,u)$, denote the total divergence

$\nabla_{x}F(x,u(x))=\frac{\partial}{\partial
x_{0}}(F^{1}(x,u(x)))+\ldots+\frac{\partial}{\partial x_{n}}(F^{n}(x,u(x))).$

When $\nabla_{x}F(x,u(x))$ is added to $L$, then $S$ can at most change by a
constant, according to the divergence theorem. Therefore, the Euler–Lagrange
equations (2) are not effected, i.e.

$\mathrm{EL}\left(L+\nabla_{x}F(x,u(x))\right)=\mathrm{EL}(L).$

######  Example 1 (Linear motions).

Denote $x_{0}=t$. Let $A\in\mathbb{R}^{d\times d}$ be a non-degenerate,
symmetric matrix. For $L(u,u_{t})=\frac{1}{2}u_{t}^{\top}Au_{t}$ the
Euler–Lagrange equation (2) yields $Au_{tt}=0$. If $A$ is non-degenerate, the
field theory describes linear motions in $\mathbb{R}^{d}$. Indeed, for any
$\theta\colon\mathbb{R}^{d}\to\mathbb{R}$ with nowhere degenerate Hesse matrix
$\mathrm{Hess}\theta$, the field theory described by the Lagrangian density
$L(u,u_{t})=\theta(u_{t})$ yield the Euler–Lagrange equations

$0=\mathrm{EL}(L)=\mathrm{Hess}\theta(u_{t})u_{tt}$

that is solved by affine linear motions $u_{tt}=0$ in $\mathbb{R}^{d}$.

1 illustrates that Lagrangian densities are not uniquely determined by a
system’s motions. The non-uniqueness is not only related to gauge freedom (1)
but also various different actions can yield equivalent Euler–Lagrange
equations. This observation will become important later when we will learn
field theories from data.

We proceed with introducing our main examples of field theories.

######  Example 2 (Wave equation).

For a potential $V\colon\mathbb{R}\to\mathbb{R}$ the Euler-Lagrange equations
to the Lagrangian

$L(u,u_{t},u_{x})=\frac{1}{2}(u_{t}^{2}-u_{x}^{2})-V(u)$ (3)

yield the non-linear wave equation

$u_{tt}(t,x)-u_{xx}(t,x)+\nabla V(u(t,x))=0.$ (4)

Here $\nabla V$ denotes the gradient of $V$ and
$u_{tt}=\frac{\partial^{2}u}{\partial t^{2}}$,
$u_{xx}=\frac{\partial^{2}u}{\partial x^{2}}$ partial derivatives.

######  Example 3 (Degenerate Lagrangian).

Let $J=\begin{pmatrix}0&-I_{n}\\\ I_{n}&0\end{pmatrix}$, where
$I_{n}\in\mathbb{R}^{n\times n}$ is an identity matrix. For a function
$H\colon\mathbb{R}^{2n}\to\mathbb{R}$ and the Lagrangian

$L(u,u_{t})=\frac{1}{2}(J^{-1}u)^{\top}u_{t}-H(u)$ (5)

the Euler-Lagrange equations are

$Ju_{t}=\nabla H(u).$ (6)

Equation (6) are Hamilton’s equations on $\mathbb{R}^{2n}$ to the Hamiltonian
$H$, with symplectic structure on $\mathbb{R}^{2n}$ defined by $J$.

The Lagrangian (5) is degenerate, since the matrix
$\frac{\partial^{2}L}{\partial u_{t}^{2}}=0$ is not regular. Indeed, the
Euler–Lagrange equations (6) are first order differential equations rather
than of second order, as one might expect from (2).

######  Remark 2.

The Hamiltonian system (6) in 3 for the degenerate Lagrangian $L(u,u_{t})$ of
(5) may not be confused with Hamiltonian systems induced by non-degenerate
Lagrangians $\tilde{L}(u,u_{t})$ with regular matrix
$\frac{\partial^{2}\tilde{L}}{\partial u_{t}^{2}}$. These are defined on a
different phase space: consider $z=(u,p)$. Let $\tilde{n}$ denote the
dimension of the variable $u$ and let
$\tilde{J}=\begin{pmatrix}0&-I_{\tilde{n}}\\\ I_{\tilde{n}}&0\end{pmatrix}$.
Motions defined by $\mathrm{EL}(\tilde{L})=0$ (see (2)) fulfil
$\tilde{J}z=\nabla\tilde{H}(z)$ for the Hamiltonian
$\tilde{H}(z)=p^{\top}u_{t}-\tilde{L}(u,u_{t})$. Here $u_{t}$ is expressed in
terms of $z=(u,p)$ using $p=\frac{\partial\tilde{L}}{\partial
u_{t}}(u,u_{t})$.

Using a natural extension of the described variational principle (1) to
complex valued fields, we present a variational principle for Schrödinger
equationDeriglazov (2009).

######  Example 4 (Schrödinger equation).

For a potential $V\colon\mathbb{R}\to\mathbb{R}$ the Euler-Lagrange equations
to the Lagrangian

$L(\Psi,\Psi_{t},\Psi_{x})=-\frac{\hbar}{2}(\overline{\Psi}\Psi_{t}-\Psi\overline{\Psi_{t}})-\overline{\Psi_{x}}\Psi_{x}-V(\overline{\Psi}\Psi)$
(7)

yield the non-linear Schrödinger equation

$\hbar\mathrm{i}\Psi_{t}=(-\Delta+V^{\prime}(|\Psi|^{2}))\Psi.$ (8)

Here $\Psi$ is a complex valued function defined on a space-time domain
parametrised by the variables $(t,x)$, $\overline{\Psi}$ its complex
conjugation, $\hbar$ is the Plank constant, $\mathrm{i}$ the imaginary unit,
and $\Delta=\frac{\partial^{2}}{\partial x^{2}}$ denotes the spatial
Laplacian.

Alternatively, with $\Psi=\phi+\mathrm{i}p$,
$u=\begin{pmatrix}\phi,&p\end{pmatrix}^{\top}$, $J$ as in 3, the Schrödinger
equation can be obtained from the a real variational principle with Lagrangian

$L(u,u_{t},u_{x})=\hbar(J^{-1}u)^{\top}u_{t}-\|u_{x}\|^{2}-V(\|u\|^{2})$ (9)

which yields

$\hbar Ju_{t}=(-\Delta+V^{\prime}(\|u\|^{2}))u,$ (10)

which is equivalent to (8).

###  II.2 Discrete variational principles

Let us introduce discrete actions based on various types of discrete
Lagrangian densities. These can be interpreted as approximations to exact
discrete Lagrangians (see Appendix A). While we provide a rather general setup
and expression of discrete Lagrangian densities in the following, 6 introduces
the types of densities used in our numerical experiments. For a detailed
introduction we refer to the articleMarsden  _et al._ (2001) by Marsden,
Pekarsky, Shkoller, West.

First we introduce some notation for meshes. Consider

  * •

the $n+1$-dimensional cube
$X=[x^{0}_{0},x^{N_{0}}_{0}]\times\ldots\times[x^{0}_{n},x^{N_{n}}_{n}]\subset\mathbb{R}^{n+1}$,

  * •

the mesh
$X_{d}=\\{x^{0}_{0},\ldots,x^{N_{0}}_{0}\\}\times\ldots\times\\{x^{0}_{n},\ldots,x^{N_{n}}_{n}\\}\subset
X$,

  * •

the grid’s interior points
$\mathring{X}_{d}=\\{x^{1}_{0},\ldots,x^{N_{0}-1}_{0}\\}\times\ldots\times\\{x^{1}_{n},\ldots,x^{N_{n}-1}_{n}\\}\subset
X_{d}$.

Assume for simplicity that the grid $X_{d}$ is uniform with spacing $\Delta
x_{r}=(x^{N_{r}}_{r}-x^{0}_{r})/N_{r}$ ($r=0,\ldots,n$). Let $\Delta x=\Delta
x_{0}\cdot\ldots\cdot\Delta x_{n}$ denote the discrete volume element.
Consider

  * •

the index set
$\mathcal{I}=\\{0,\ldots,N_{0}\\}\times\ldots\times\\{0,\ldots,N_{n}\\}$,

  * •

the multi-index $l=(l_{0},\ldots,l_{n})\in\mathcal{I}$,

  * •

the small cube
$X^{l}=[x_{0}^{l_{0}},x_{0}^{l_{0}+1}]\times\ldots\times[x_{n}^{l_{n}},x_{n}^{l_{n}+1}]$
to $l\in{\mathcal{I}}$ with $l_{r}<N_{r}$ for all $0\leq r\leq n$,

  * •

and the set of vertices
$X_{{\mathrm{vrtx}}}^{l}=\\{x_{0}^{l_{0}},x_{0}^{l_{0}+1}\\}\times\ldots\times\\{x_{n}^{l_{n}},x_{n}^{l_{n}+1}\\}$
of a small cube $X^{l}$. Each $X_{{\mathrm{vrtx}}}^{l}$ has $2^{n+1}$
elements.

A discrete action functional $S_{d}$ now assigns real values to functions
$U\colon X_{d}\to\mathbb{R}$ defined on the grid $X_{d}$: a discrete action
functional is given as $S_{d}\colon(\mathbb{R}^{d})^{X_{d}}\to\mathbb{R}$ with

$S_{d}(U)=\sum_{l\in\mathcal{I},l_{r}<N_{r}}L_{d}(X^{l}_{{\mathrm{vrtx}}},U(X^{l}_{{\mathrm{vrtx}}}))\Delta
x,$ (11)

where
$L_{d}\colon(\mathbb{R}^{n+1})^{2^{n+1}}\times(\mathbb{R}^{d})^{2^{n+1}}\to\mathbb{R}$
can be evaluated from the values of $U$ on each vertex set
$X^{l}_{\mathrm{vrtx}}$. $L_{d}$ is called discrete Lagrangian density. In
many cases $L_{d}$ is autonomous, i.e. does not explicitly depend on
$X^{l}_{\mathrm{vrtx}}$ and is, thus, of the form
$L_{d}(U(X^{l}_{{\mathrm{vrtx}}}))$. The discrete action $S_{d}$ of (11) is a
discretised version of (1).

The discrete Euler-Lagrange equations are obtained by discrete variations
fixing the boundary as $\nabla_{U(\mathring{X}_{d})}S_{d}(U)=0$. This gives
for any interior vortex $v\in\mathring{X_{d}}$ the discrete Euler–Lagrange
equation

$\mathrm{DEL}(L_{d})_{v}:=\nabla_{U(v)}\sum_{l\in\mathcal{I}_{v}}L_{d}(X^{l}_{{\mathrm{vrtx}}},U(X^{l}_{{\mathrm{vrtx}}}))=0,$
(12)

where $\mathcal{I}_{v}$ contains all indices $l$ of cubes that contain $v$,
i.e. $v\in X^{l}$. The derivative $\nabla_{U(v)}$ describes the derivative
with respect to the $u$-variable at the vortex $v$.

######  Example 5 (Linear motions).

Denote $x_{0}=t$. Let $A\in\mathbb{R}^{d\times d}$ be a symmetric, non-
degenerate matrix. Consider a uniform mesh $X_{d}=\\{t^{0},\ldots,t^{N}\\}$ on
$[t^{0},t^{N}]\subset\mathbb{R}$. We identify $(\mathbb{R}^{d})^{X_{d}}$ with
the space $(\mathbb{R}^{d})^{N+1}$. The discrete Lagrangian density

$L_{d}(u^{i},u^{i+1})=\frac{1}{2}(u^{i+1}-u^{i})^{\top}A(u^{i+1}-u^{i})$

yields discrete Euler–Lagrange equations

$\displaystyle 0=\mathrm{DEL}(L_{d})_{u^{i}}$
$\displaystyle=\nabla_{u^{i}}\left(L_{d}(u^{i-1},u^{i})+L_{d}(u^{i},u^{i+1})\right)$
$\displaystyle=A(2u^{i}-u^{i-1}-u^{i+1}).$

These are uniquely solved by discrete linear dynamics

$u^{i+1}=u^{i}+(u^{i}-u^{i-1}).$ (13)

Indeed, in analogy to 1, for any $\theta\colon\mathbb{R}^{d}\to\mathbb{R}$
with nowhere degenerate Hesse matrix the discrete Lagrangian

$L_{d}(u^{i},u^{i+1})=\theta(u^{i+1}-u^{i})$

yields discrete Euler–Lagrange equations

$0=\mathrm{DEL}(L_{d})_{u^{i}}=\nabla\theta(u^{i}-u^{i-1})-\nabla\theta(u^{i+1}-u^{i})$

to which (13) is the locally unique solution.

5 demonstrates that discrete Lagrangian densities are not uniquely determined
by a system’s motion. Depending on the choice of $L_{d}$, the discrete
Euler–Lagrange equations can be easier or more difficult to solve numerically.
This will need to be considered when learning a model of $L_{d}$ from data.
The following example introduces the two main stencils used in our numerical
experiments.

Figure 2: Visualisation of the 7 and 9 point stencil of 6. Colours indicate
that points are part of the same summand in (14) or (15). Black vertices are
present in several summands.

######  Example 6 (3 and 4 point discrete Lagrangian).

We denote $x_{0}=t$, $x_{1}=x$ and identify
$\mathbb{R}^{X_{d}}\cong\mathbb{R}^{(N_{0}+1)\times(N_{1}+1)}$ by the space of
$(N_{0}+1)\times(N_{1}+1)$ matrices with elements $U=(u^{i}_{j})^{0\leq i\leq
N_{0}}_{0\leq j\leq N_{1}}$. For a discrete Lagrangian of the form
$L_{d}(u^{i}_{j},u^{i+1}_{j},u^{i}_{j+1})$ (3 input arguments, autonomous) the
discrete Euler-Lagrange equations (12) yield the 7 point stencil

$\begin{split}\mathrm{DEL}(L_{d})_{j}^{i}=\frac{\partial}{\partial
u^{i}_{j}}\big{(}&\phantom{+}L_{d}(u^{i}_{j},u^{i+1}_{j},u^{i}_{j+1})\\\
&+L_{d}(u^{i-1}_{j},u^{i}_{j},u^{i-1}_{j+1})\\\
&+L_{d}(u^{i}_{j-1},u^{i+1}_{j-1},u^{i}_{j})\big{)}=0,\end{split}$ (14)

denoting $\mathrm{DEL}(L_{d})_{j}^{i}=\mathrm{DEL}(L_{d})_{u^{i}_{j}}$. If
$L_{d}$ is of the form
$L_{d}(u^{i}_{j},u^{i+1}_{j},u^{i}_{j+1},u^{i+1}_{j+1})$ (4 input arguments,
autonomous), then the discrete Euler-Lagrange equations yield the 9 point
stencil

$\begin{split}\mathrm{DEL}(L_{d})_{j}^{i}=\frac{\partial}{\partial
u^{i}_{j}}\big{(}&\phantom{+}L_{d}(u^{i}_{j},u^{i+1}_{j},u^{i}_{j+1},u^{i+1}_{j+1})\\\
&+L_{d}(u^{i-1}_{j},u^{i}_{j},u^{i-1}_{j+1},u^{i}_{j+1})\\\
&+L_{d}(u^{i}_{j-1},u^{i+1}_{j-1},u^{i}_{j},u^{i+1}_{j})\\\
&+L_{d}(u^{i-1}_{j-1},u^{i}_{j-1},u^{i-1}_{j},u^{i}_{j})\big{)}=0.\end{split}$
(15)

The stencils are visualised in Fig. 2

######  Remark 3 (Gauge freedom).

In analogy to 1, discrete Lagrangians are not uniquely determined by a
system’s discrete motions. Indeed, two 3-point Lagrangians $L_{d}(a,b,c)$ and
$\tilde{L}_{d}(a,b,c)$ as in 6 yield the same discrete Euler–Lagrange
equations (14) if they differ by a discrete divergence, i.e. if
$L_{d}(a,b,c)-\tilde{L}_{d}(a,b,c)$ is of the form

$\chi_{1}(a)-\chi_{1}(b)+\chi_{2}(a)-\chi_{2}(c)+\chi_{3}(b)-\chi_{3}(c)$ (16)

for a differentiable function $\chi=(\chi_{1},\chi_{2},\chi_{3})\colon
X\to\mathbb{R}^{3}$. Analogous observations apply to other types of discrete
Lagrangians.

We proceed by introducing the two main discrete field theories used in the
numerical experiments.

######  Example 7 (Discrete 2d wave equation).

We denote $x_{0}=t$, $x_{1}=x$ and identify
$\mathbb{R}^{X_{d}}\cong\mathbb{R}^{(N_{0}+1)\times(N_{1}+1)}$ by the space of
$(N_{0}+1)\times(N_{1}+1)$ matrices with elements $U=(u^{i}_{j})^{0\leq i\leq
N_{0}}_{0\leq j\leq N_{1}}$. Let $L$ be the Lagrangian of 2. Consider the
autonomous discrete Lagrangian density

$\begin{split}L_{d}(u^{i}_{j},u^{i+1}_{j},u^{i}_{j+1})=L\left(u^{i}_{j},\frac{u^{i+1}_{j}-u^{i}_{j}}{\Delta
t},\frac{u^{i}_{j+1}-u^{i}_{j}}{\Delta x}\right)\\\
=\frac{1}{2}\left(\frac{u^{i+1}_{j}-u^{i}_{j}}{\Delta
t}\right)^{2}-\frac{1}{2}\left(\frac{u^{i}_{j+1}-u^{i}_{j}}{\Delta
x}\right)^{2}-V(u^{i}_{j}).\end{split}$ (17)

The discrete Euler–Lagrange equation (14) is given as

$\frac{u^{i-1}_{j}-2u^{i}_{j}+u^{i+1}_{j}}{\Delta
t^{2}}-\frac{u^{i}_{j-1}-2u^{i}_{j}+u^{i}_{j+1}}{\Delta x^{2}}+\nabla
V(u^{i}_{j})=0.$ (18)

It can be interpreted as a 2nd order discretisation of 2 by finite
differences. While for general $L_{d}$ the discrete Euler–Lagrange equation
(14) constitutes a 7-point stencil, on $L_{d}$ of (17) it simplifies to a
5-point stencil. The stencil can be used conveniently to compute forward
propagations of solutions as $u^{i+1}_{j}$ can be computed in terms of
$u^{i}_{j}$, $u^{i-1}_{j}$, $u^{i}_{j+1}$, which all belong to previous time-
steps.

The following example uses a straight forward extension of the discrete
variational principle to complex valued fields.

######  Example 8 (Discrete Schrödinger equation).

Let $L$ be the Lagrangian (7) of 4. Consider

$L_{d}(\Psi^{i-1}_{j-1},\Psi^{i}_{j-1},\Psi^{i-1}_{j},\Psi^{i}_{j})=L(m^{i}_{j}\Psi,\Delta
t^{i}_{j}\Psi,\Delta x^{i}_{j}\Psi)$ (19)

with

$\displaystyle m^{i}_{j}\Psi$
$\displaystyle=\frac{\Psi^{i-1}_{j-1}+\Psi^{i-1}_{j}+\Psi^{i}_{j}+\Psi^{i}_{j-1}}{4}$
(20) $\displaystyle\Delta t^{i}_{j}\Psi$
$\displaystyle=\frac{(\Psi^{i}_{j-1}-\Psi^{i-1}_{j-1})+(\Psi^{i}_{j}-\Psi^{i-1}_{j})}{2\Delta
t}$ (21) $\displaystyle\Delta x^{i}_{j}\Psi$
$\displaystyle=\frac{(\Psi^{i-1}_{j}-\Psi^{i-1}_{j-1})+(\Psi^{i}_{j}-\Psi^{i}_{j-1})}{2\Delta
x}$ (22)

that approximate the field $\Psi$ and its derivative in $t$ and $x$ in the
centre of a square $X^{(i,j)}$ of the gird. Further, consider the following
approximations of derivatives at the node $(i,j)$

$\displaystyle{\mathrm{D}_{t}}^{i}_{j}\Psi$
$\displaystyle=\frac{(\Psi^{i+1}_{j-1}-\Psi^{i-1}_{j-1})+2(\Psi^{i+1}_{j}-\Psi^{i-1}_{j})+(\Psi^{i+1}_{j+1}-\Psi^{i-1}_{j+1})}{8\Delta
t}$ $\displaystyle{\mathrm{D}_{x}^{2}}^{i}_{j}\Psi$
$\displaystyle=\frac{1}{4\Delta
x^{2}}\Big{(}(\Psi^{i-1}_{j-1}-2\Psi^{i-1}_{j}+\Psi^{i-1}_{j+1})$
$\displaystyle\phantom{=\frac{1}{4\Delta
x^{2}}\Big{(}}+(\Psi^{i}_{j-1}-2\Psi^{i}_{j}+\Psi^{i}_{j+1})$
$\displaystyle\phantom{=\frac{1}{4\Delta
x^{2}}\Big{(}}+(\Psi^{i+1}_{j-1}-2\Psi^{i+1}_{j}+\Psi^{i+1}_{j+1})\Big{)}.$

As discrete Euler-Lagrange equations (15) we obtain

$\begin{split}\mathrm{i}\hbar{\mathrm{D}_{t}}^{i}_{j}\Psi&=-{\mathrm{D}_{x}^{2}}^{i}_{j}\Psi\\\
&+\frac{1}{4}\Big{(}V^{\prime}(m^{i-1}_{j-1}\Psi)m^{i-1}_{j-1}\Psi+V^{\prime}(m^{i}_{j-1}\Psi)m^{i}_{j-1}\Psi\\\
&+\phantom{\frac{1}{4}\Big{(}}V^{\prime}(m^{i-1}_{j}\Psi)m^{i-1}_{j}\Psi+V^{\prime}(m^{i}_{j}\Psi)m^{i}_{j}\Psi\Big{)},\end{split}$
(23)

which is a 9 point stencil.

######  Remark 4 (Computations).

In 6 assume that boundary conditions $(u_{0}^{i})^{0\leq i\leq N_{0}}$,
$(u_{N_{1}}^{i})^{0\leq i\leq N_{0}}$, initial data $(u_{j}^{0})_{0<j<N_{1}}$,
and initial velocities $(\dot{u}_{j}^{0})_{0<j<N_{1}}$ are given. Let

$L_{\Delta x}^{\Delta
t}(U^{i},U^{i+1})=\sum_{j=0}^{N_{1}-1}L_{d}(u^{i}_{j},u_{j}^{i+1},u^{i}_{j+1})$
(24)

or

$L_{\Delta x}^{\Delta
t}(U^{i},U^{i+1})=\sum_{j=0}^{N_{1}-1}L_{d}(u^{i}_{j},u_{j}^{i+1},u^{i}_{j+1},u^{i+1}_{j+1}),$
(25)

respectively. Here $U^{i}=(u^{i}_{j})_{0<j<N_{1}}$ refers to values at
interior spatial grid points. Let $L_{\Delta x}(U,\dot{U})=L_{\Delta
x}^{\Delta t}(U-\Delta t/2\dot{U},U+\Delta t/2\dot{U})$. Solutions to the
discrete field theory can be computed as follows: $U^{1}$ is computed as the
solution to

$\nabla_{\dot{U}}L_{\Delta x}(U^{0},\dot{U}^{0})=-\nabla_{U^{1}}L_{\Delta
x}^{\Delta t}(U^{0},U^{1}).$ (26)

Subsequent $U^{i}$ are computed as solutions to the DEL for $L_{\Delta
x}^{\Delta t}$

$\nabla_{U^{i}}\left(L_{\Delta x}^{\Delta t}(U^{i-1},U^{i})+L_{\Delta
x}^{\Delta t}(U^{i},U^{i+1})\right)=0.$ (27)

This follows classical variational integration theoryMarsden and West (2001).
We refer to the literature on low-rank approximations and integration
schemesKieri, Lubich, and Walach (2016); Lubich, Vandereycken, and Walach
(2018); Walach (2019); Schrammer (2022) in case fine meshes or higher spacial
dimensions are considered and the systems (26), (27) become prohibitively
high-dimensional.

In the special case that $L_{\Delta x}(U,\dot{U})$ is linear in the velocities
$\dot{U}$, (26) is independent of $\dot{U}^{0}$ such that no initial velocity
profile $\dot{U}^{0}$ is needed to obtain $U^{1}$ via (26). In particular, if
$L_{\Delta x}(U,\dot{U})$ is of the form (5), then the computation is
numerically stableRowley and Marsden (2002). Indeed, independence of
$\dot{U}^{0}$ of the numerical scheme is expected since the EL equations (6)
are 1st order ordinary differential equations. While numerical stability of
the scheme is guaranteed for $L_{\Delta x}$ of the form (5), for more general
Lagrangians linear in velocities other techniques are required for numerically
stable computationsKraus (2017); Ellison  _et al._ (2018); Burby, Finn, and
Ellison (2022).

######  Remark 5.

Depending on boundary conditions, the solution procedure of 4 can often be
drastically simplified by using the local stencils (14), (15) to propagate
solutions forward in time, as we made explicit in 7. Use of the local stencils
avoids solving the potentially high dimensional system (27).

######  Remark 6 (Numerical stability of discrete SE).

Solutions to (23) can be computed as explained in 4. Notice that when the
discrete Lagrangian $L_{d}$ of (19) is expressed in the variable $(u,\phi)$,
where $\Psi=u+\mathrm{i}\phi$, then $L_{\Delta x}^{\Delta t}$ of 4 is of the
form (5) of 3. Therefore, the computation can be initialised from initial data
$(\Psi_{j}^{0})_{0\leq j\leq N_{1}}$ and boundary conditions
$(\Psi_{0}^{i})^{0<i\leq N_{0}}$, $(\Psi_{N_{1}}^{i})^{0<i\leq N_{0}}$ and is
numerically stableRowley and Marsden (2002). No velocity data is required for
the initialisation.

##  III Preparation of regularisation strategy

In the following, we will train a neural network model of a discrete
Lagrangian density $L_{d}$ such that $L_{d}$ defines a discrete field theory
that is consistent with observations. Discrete Lagrangian densities are not
uniquely determined by the system’s motions as illustrated in 5. (Also see 3.)
The secondary goal is to identify a model of $L_{d}$ that is not only
consistent with the true dynamics but is also computationally efficient when
numerical methods are used to solve the discrete Euler–Lagrange equations.

To prepare our design of regularisers to be employed in the training phase, we
analyse the convergence properties of Newton interations when used to compute
solutions to discrete Euler–Lagrange equations for the 3 and 4 point discrete
Lagrangians of 6. The following Section refers to the DEL for 3-point
Lagrangians (14).

######  Proposition 1.

Let $u^{i}_{j}$, $u^{i+1}_{j}$, $u^{i}_{j+1}$, $u^{i-1}_{j}$, $u^{i-1}_{j+1}$,
$u^{i}_{j-1}$, $u^{i+1}_{j-1}$ such that (14) holds. Let
$\mathcal{O}\subset\mathbb{R}^{d}$ be a convex, neighbourhood of
$u^{\ast}=u^{i+1}_{j}$, $\|\cdot\|$ a norm of $\mathbb{R}^{d}$ inducing an
operator norm on $\mathbb{R}^{d\times d}$. Define
$p(u):=\frac{\partial^{2}L_{d}}{\partial u^{i}_{j}\partial
u}(u^{i}_{j},u,u^{i}_{j+1})$ and let $\theta$ and $\overline{\theta}$ be
Lipschitz constants on $\mathcal{O}$ for $p$ and for $\mathrm{inv}\circ p$,
respectively, where $\mathrm{inv}$ denotes matrix inversion. Let

$\rho^{\ast}:=\left\|\mathrm{inv}(p(u^{\ast}))\right\|=\left\|\left(\frac{\partial^{2}L_{d}}{\partial
u^{i}_{j}\partial
u^{\ast}}(u^{i}_{j},u^{\ast},u^{i}_{j+1})\right)^{-1}\right\|$ (28)

and let $f(u^{(n)})$ denote the left hand side of (14) with $u_{j}^{i+1}$
replaced by $u^{(n)}$. If
$\|{u}^{(0)}-{u}^{\ast}\|\leq\min\left(\frac{\rho^{\ast}}{\overline{\theta}},\frac{1}{2\theta\rho^{\ast}}\right)$
for ${u}^{(0)}\in\mathcal{O}$, then the Newton Iterations
${u}^{(n+1)}:={u}^{(n)}-\mathrm{inv}(p(u^{(n)}))f(u^{(n)})$ converge
quadratically against ${u}^{\ast}$, i.e.

$\|{u}^{(n+1)}-{u}^{\ast}\|\leq\rho^{\ast}\theta\|{u}^{(n)}-{u}^{\ast}\|^{2}.$
(29)

A proof can be found in Appendix B.

######  Remark 7.

1 is a restatement of Proposition 1 of our conference contributionOffen and
Ober-Blöbaum (2023), where it is stated without proof.

1 clarifies numerical convergence properties when computing solutions on a
grid with the stencil visualised to the left of Fig. 2. It applies when values
$u_{0}^{i}$ at the left boundary and initial values $u_{j}^{1}$, $u_{j}^{2}$
are known for all indices $i,j$ such that remaining values can be computed
subsequently. However, sometimes it is required to compute all spatial data at
a time step at once due to boundary conditions, particularities of the
considered stencils or numerical stability issues. This is the setting of the
computations mentioned in 4. To obtain a field theory optimised for such
computations, we present the following Section.

######  Proposition 2.

Let $U^{i-1}$, $U^{i}$, $U^{i+1}\in(\mathbb{R}^{d})^{(N_{1}-1)}$ such that
(27) holds. Let $O\subset(\mathbb{R}^{d})^{(N_{1}-1)}$ be a convex,
neighbourhood of $u^{\ast}=U^{i+1}$, $\|\cdot\|$ a norm on
$(\mathbb{R}^{d})^{(N_{1}-1)}$ inducing an operator norm on
$(\mathbb{R}^{d})^{(N_{1}-1)\times(N_{1}-1)}$. Define
$p(U):=\frac{\partial^{2}L_{\Delta x}^{\Delta t}}{\partial U^{i}\partial
U}(U^{i},U)$ and let $\theta$ and $\overline{\theta}$ be Lipschitz constants
on $\mathcal{O}$ for $p$ and for $\mathrm{inv}\circ p$, respectively, where
$\mathrm{inv}$ denotes matrix inversion. Let

$\rho^{\ast}:=\left\|\mathrm{inv}(p(U^{\ast}))\right\|=\left\|\left(\frac{\partial^{2}L_{\Delta
x}^{\Delta t}}{\partial U^{i}\partial
U^{\ast}}(U^{i},U^{\ast})\right)^{-1}\right\|$ (30)

and let $f(U^{(n)})$ denote the left hand side of (27) with $U^{i+1}$ replaced
by $U^{(n)}$. If
$\|{U}^{(0)}-{U}^{\ast}\|\leq\min\left(\frac{\rho^{\ast}}{\overline{\theta}},\frac{1}{2\theta\rho^{\ast}}\right)$
for ${u}^{(0)}\in\mathcal{O}$, then the Newton Iterations
${U}^{(n+1)}:={U}^{(n)}-\mathrm{inv}(p(U^{(n)}))f(U^{(n)})$ converge
quadratically against ${U}^{\ast}$, i.e.

$\|{U}^{(n+1)}-{U}^{\ast}\|\leq\rho^{\ast}\theta\|{U}^{(n)}-{U}^{\ast}\|^{2}.$
(31)

###### Proof.

The proof is analogous to the proof of 1, which can be found in Appendix B. ∎

######  Remark 8.

The matrix $\Lambda=\frac{\partial^{2}L_{\Delta x}^{\Delta t}}{\partial
U^{i}\partial U^{i+1}}(U^{i},U^{i+1})$ of (30) has the following block-tri-
diagonal structure

$\Lambda=\begin{pmatrix}A_{1}&B_{1}\\\ C_{2}&A_{2}&B_{2}\\\
&\ddots&\ddots&\ddots\\\ &&C_{N_{1}-2}&A_{N_{1}-2}&B_{N_{1}-2}\\\
&&&C_{N_{1}-1}&A_{N_{1}-1}\\\
\end{pmatrix}\in(\mathbb{R}^{d})^{(N_{1}-1)\times(N_{1}-1)},$ (32)

where

$\begin{split}A_{j}=\frac{\partial^{2}}{\partial u^{i}_{j}\partial
u^{i+1}_{j}}\big{(}&\phantom{+}L_{d}(u^{i}_{j},u^{i+1}_{j},u^{i}_{j+1},u^{i+1}_{j+1})\\\
&+L_{d}(u^{i}_{j-1},u^{i+1}_{j-1},u^{i}_{j},u^{i+1}_{j})\big{)}\\\
B_{j}=\frac{\partial^{2}}{\partial u^{i}_{j}\partial
u^{i+1}_{j+1}}\big{(}&\phantom{+}L_{d}(u^{i}_{j},u^{i+1}_{j},u^{i}_{j+1},u^{i+1}_{j+1})\big{)}\\\
C_{j}=\frac{\partial^{2}}{\partial u^{i}_{j}\partial
u^{i+1}_{j-1}}\big{(}&\phantom{+}L_{d}(u^{i}_{j-1},u^{i+1}_{j-1},u^{i}_{j},u^{i+1}_{j})\big{)}.\end{split}$
(33)

If periodic boundary conditions rather then Dirichlet boundary conditions are
considered then the spatial mesh has $N_{1}$ interior mesh points and

$\Lambda=\begin{pmatrix}A_{0}&B_{0}&&&C_{0}\\\ C_{1}&A_{1}&B_{1}\\\
&\ddots&\ddots&\ddots\\\ &&C_{N_{1}-2}&A_{N_{1}-2}&B_{N_{1}-2}\\\
B_{N_{1}-1}&&&C_{N_{1}-1}&A_{N_{1}-1}\\\
\end{pmatrix}\in(\mathbb{R}^{d})^{N_{1}\times N_{1}}$ (34)

with $A_{j}$, $B_{j}$, $C_{j}$ as in (33)

######  Remark 9 (Numerical computation of $\|\Lambda^{-1}\|$).

Let $\|\Lambda\|$ denote the spectral norm of $\Lambda$. The value
$\rho^{\ast}=\|\Lambda^{-1}\|^{-1}$ is the smallest singular value $\sigma$ of
$\Lambda$. Its square $\sigma^{2}$ can be approximated using inverse vector
iterationsDeuflhard and Hohmann (2003a) applied to $M=\Lambda^{\top}\Lambda$:
given a random start vector $v^{\prime}$, iterate $v\leftarrow v^{\prime}$,
$v\leftarrow M\backslash v$, $v\leftarrow v/\|v\|$,
$\sigma^{2}\leftarrow(v^{\top}v^{\prime})^{-1}$. Here $M\backslash v$ denotes
the solution of the linear equation $Mx=v$. For this a Cholesky factorisation
of $M$ can be computed and used throughout the iteration. While in our the
case sparse structure of $\Lambda$ can be exploited, we refer to the
literatureDeuflhard and Hohmann (2003a) for methods that avoid the explicit
computation of the matrix product $\Lambda^{\top}\Lambda$ in the computation
of singular values. Further, more elaborate tensor representationsRöhrig-
Zöllner, Thies, and Basermann (2022) can be used, especially when generalising
to higher-dimensional space-time lattices. Approximations to
$\sigma^{2}=(\rho^{\ast})^{2}$ will be used in the machine learning framework
introduced in the following sections.

##  IV Loss function and regulariser for training discrete Lagrangian
densities

Assume we are given observations $U^{(k)}\colon X_{d}\to\mathbb{R}^{d}$
$k=1,\ldots,K$ of solutions to an unknown discrete field theory and the goal
is to identify a model of the discrete field theory. We pick a type of a
discrete Lagrangian density $L_{d}$ such as the 3 or 4 points discrete
Lagrangians introduced in 6 or any discretisation of (55) by a linear
quadrature formula. The discrete Lagrangian $L_{d}$ is modelled as a neural
network whose parameters are optimised such that

  1. 1.

the discrete Euler–Lagrange equations (12) are consistent with the
observations

  2. 2.

and such that computations based on (12) are numerically well conditioned.

For this, we consider a loss function
$\ell=\ell_{\mathrm{data}}+w\ell_{\mathrm{reg}}$ consisting of the data
consistency term $\ell_{\mathrm{data}}$ and the regulariser
$\ell_{\mathrm{reg}}$ and a weight $w\in\mathbb{R}_{+}$. These are specified
in the following.

###  IV.1 Data consistency

We introduce the data consistency term

$\ell_{\mathrm{data}}=\sum_{k=1}^{K}\sum_{v}|\mathrm{DEL}(L_{d})_{v}(U^{(k)})|^{2}.$
(35)

In the inner sum, the sum is taken over all vertices for which the summand is
defined. The term $\ell_{\mathrm{data}}$ measures how consistent the model
$L_{d}$ is with the training data.

######  Remark 10 (Partial observations and batch learning).

Since the computation of a summand $|\mathrm{DEL}(L_{d})_{v}(U^{(k)})|^{2}$ in
(35) does not involve the full observation $U^{(k)}(X_{d})$ but only the
observation on a stencil such as the ones illustrated in Fig. 2 for 3- and
4-point discrete Lagrangian densities $L_{d}$, partial observations can be
incorporated in $\ell_{\mathrm{data}}$. This is especially useful when batch
learningMarsland (2014) is applied.

######  Remark 11 (Low dimensionality).

Although solutions to the discrete field theory can require solving high
dimensional systems such as $|\mathrm{DEL}(L_{\Delta x}^{\Delta
t})^{i}(U)|^{2}=0$ in (27), in the training phase to compute
$\ell_{\mathrm{data}}$ each summand only requires evaluating the low-
dimensional Euler–Lagrange equations (24).

######  Remark 12 (Locality).

In the setting of the article, we assume that there exists a continuous field
theory underlying our observational data $U^{(k)}$ on a mesh $X_{d}$
(spacetime lattice). The underlying continuous theory can be described by a
partial differential equation (2) and whether a function $u\colon
X\to\mathbb{R}$ fulfils the differential equation can be checked point-wise,
i.e. it is a local condition. This locality is reflected in our discrete model
by the locality of the stencil defining the term
$\mathrm{DEL}(L_{d})_{v}(U^{(k)})$. In this sense, the data-driven model will
be local in space and time by construction. In contrast, if in the setting of
a space-time lattice $X=[t^{0},t^{N}]\times[x^{0},x^{N}]$ we were to learn
$L_{\Delta x}^{\Delta t}$ (see (25)) rather than the discrete Lagrangian
density $L_{d}$, potentially after model order reduction in the spatial
directions, then locality would only be guaranteed in the time-direction.

###  IV.2 Regulariser

Since any constant function $L_{d}$ constitutes a minimiser of
$\ell_{\mathrm{data}}$, a carefully chosen regulariser $\ell_{\mathrm{reg}}$
is required to obtain a regular discrete field theory. Indeed, as seen in 5
much freedom in the choice of $L_{d}$ can be expected. In view of Section 1
and 2, we consider regularisers which aim to minimise $\rho^{\ast}$ of (29) or
(31), respectively. If successfully optimised, we can expect solvability and
rapid convergence of Newton Iterations when computing solutions to the field
theory, at least for initial data close to training data.

A regulariser $\ell_{\mathrm{reg}}$ for the 3-point discrete Lagrangian of 6
is given by

$\ell_{\mathrm{reg}}=\alpha\sum_{k=1}^{K}\sum_{i,j}\left\|\left(\frac{\partial^{2}L_{d}}{\partial
u^{i}_{j}\partial
u^{i+1}_{j}}({u^{i}_{j}}^{(k)},{u^{i+1}_{j}}^{(k)},{u^{i}_{j+1}}^{(k)})\right)^{-1}\right\|^{2}$
(36)

with normalisation $\alpha=\frac{1}{KN_{0}N_{1}}$ by the number of summands.
In the summation in (36) the indices $i,j$ run over all interior mesh points.
The expression is motivated by 1: each summand corresponds to
$(\rho^{\ast})^{2}$ of the Section which is related to the convergence speed
of Newton-Iterations. We obtain the bound
$\rho^{\ast}<\sqrt{\ell_{\mathrm{reg}}}$. The resulting field theory will be
optimised for propagating solutions forward in space and time using the
7-point stencil Eq. 14 visualised to the left of Fig. 2.

To obtain a field theory based on the 3- or 4-point discrete Lagrangian of 6
which is optimised for forward propagation of solutions as described in 4, we
consider the regularisation

$\ell_{\mathrm{reg}}=\frac{1}{KN_{0}}\sum_{k=1}^{K}\sum_{i=1}^{N_{0}}\left\|{\Lambda_{i}^{(k)}}^{-1}\right\|^{2},$
(37)

with

$\Lambda^{(k)}_{i}=\frac{\partial^{2}L_{\Delta x}^{\Delta t}}{\partial
U^{i}\partial U^{i+1}}({U^{i}}^{(k)},{U^{i+1}}^{(k)})$ (38)

The expression is based on 2. Here
${U^{i}}^{(k)}=({u^{i}_{j}}^{(k)})_{0<j<N_{1}}$. See 8 and above for a
description of $\Lambda^{(k)}$. In practise, we found that the tamed version

$\ell_{\mathrm{reg}}=\frac{1}{KN_{0}}\sum_{k=1}^{K}\sum_{i=1}^{N_{0}}\mathrm{relu}\left(-10\left\|{\Lambda^{(k)}}^{-1}\right\|^{-2}+1\right)\leq
1$ (39)

is more suitable for training. Here $\mathrm{relu}(x)=\max(0,x)$ is the
rectified linear unit function. For numerical computation of
$\left\|{\Lambda^{(k)}}^{-1}\right\|^{-1}$ see 9. The regularisers (37) and
(39) both get large when the singular values of $\Lambda^{(k)}$ become small
on the training data. However, (39) is tamed in the sense that the summands
take values in $[0,1)$ and are zero when the smallest singular values of
$\Lambda^{(k)}$ is bigger than or equal to $0.1$. On the training data set,
$\rho^{\ast}$ of 2 is bounded by
$\rho^{\ast}\leq\sqrt{10(1-\ell_{\mathrm{reg}})^{-1}}$. Figure 3 shows the
behaviour of the summands of various regularisers in terms of the square of
the smallest singular value $\sigma^{2}$.

Figure 3: Behaviour of summands in regularisers $\ell_{\mathrm{reg}}$ of (36),
(37), (39), and $\ell_{\mathrm{reg}}^{\mathrm{illustration}}$.

Figure 4: Illustration of degeneracy of roots when learning without (left) or
with (right) regularising terms. The regulariser makes sure that the roots are
non-degenerate and numerical rootfinding is well-conditioned. For comparison,
a posterior distribution for a Gaussian Processes for the same observations is
computed without regularisation (right)  Figure 5: As comparison to Fig. 4
(left), a posterior distribution for a Gaussian Processes for the same
observations is computed without regularisation. The problem of degenerate
roots of the mean is less pronounced but does occur for sparse training data.

For a better intuition of the action of the considered regularisers, Fig. 4
shows the action of a regulariser in a simplified learning task: a neural
network model of a function $f\colon\mathbb{R}\to\mathbb{R}$ is fitted to
observations $f(x)=0$ for $x\in R=\\{-8,-5,1.2,3\\}$. While a model trained
with the mean squared error
$\ell_{\mathrm{data}}^{\mathrm{illustration}}=\sum_{x\in R}f(x)^{2}$ is almost
flat after $1100$ optimisation steps with adamKingma and Ba (2015) (left of
Fig. 4), the presence of the regulariser
$\ell_{\mathrm{reg}}^{\mathrm{illustration}}=\sum_{x\in
R}(0.01+f^{\prime}(x)^{2})^{-1}$ drives the model of $f$ to a representation
which has nondegenerate roots at $R$ that are easy to find using Newton
Iterations (right of Fig. 4). Indeed, the loss of the regulariser
$\ell_{\mathrm{reg}}^{\mathrm{illustration}}<0.855$ guarantees that
$f^{\prime}(x)>3.38$ for $x\in R$ and thus the regularity of the roots problem
for $f$ close to $R$ is guaranteed.

######  Remark 13 (Comparison with Gaussian Processes).

When Gaussian Processes are used instead of neural networks, the problem of
degenerate roots are not as pronounced: Fig. 5 shows the posterior
distribution and its mean of a Gaussian Process for observations $f(x)=0$,
$x\in R$ (constant prior $m_{0}(x)=0.5$, squared exponential kernel with unit
length scale and unit standard deviation). The roots of the mean function at
$1.2$ and $3$ are non-degenerate, while roots at $-8$ and $-5$ appear
tangential. This behaviour is quite distinct from that of neural networks Fig.
4 (left). We conclude that each learning approaches requires adapted
regularisation. We have developed regularisation strategies for Lagrangians
modelled as Gaussian Processes in previous workOber-Blöbaum and Offen (2023).

##  V Variational principles for travelling waves in discrete and continuous
theories

This section discusses the notion of travelling waves in discrete field
theories. As reference solutions for our numerical experiments in Section VI,
we report analytic expressions for travelling waves of the discrete wave
equation and discrete Schrödinger equation. Further, a central motivation for
learning (discrete) field theories by learning a (discrete) Lagrangian density
is the presence of variational structure of the learned model. We invoke
Palais’ principle of symmetric criticalityPalais (1979) which implies that
variational structure is passed on to the governing equations of motions of
travelling waves if the discrete Lagrangian is autonomous. In the broader
context of learning field theories and trying to predict highly symmetric
solutions with the trained model, this is used as a motivation for learning
discrete Lagrangian densities with the correct invariances.

Highly symmetric solutions $u$ of a continuous field theory, such as
travelling waves, are of special interest due to their simple structure. When
the same dynamical system is described by a discrete field theory, however,
the values $u(X_{d})$ of $u$ over the mesh $X_{d}$ might not have symmetry
properties unless the symmetries of $u$ happen to align with the symmetries of
the mesh $X_{d}$. For instance, let $X=\mathbb{R}\times[0,b]/\sim$. Here
$[0,b]/\sim$ denotes the identification of endpoints $0\sim b$ of an interval
$[0,b]$ (periodic boundary conditions). Functions $u\colon X\to\mathbb{R}$
with $u(t,x)=u(t+s,x+cs)$ for all $s\in\mathbb{R}$ are of the form
$u(t,x)=f(x-ct)$ for $f\colon[0,b]/\sim\to\mathbb{R}$ (travelling wave with
wave speed $c\in\mathbb{R}$). Consider a rectangular mesh $X_{d}\subset X$
with mesh widths $\Delta t$, $\Delta x>0$. If and only if $\frac{c\Delta
t}{\Delta x}\in\mathbb{Q}$, the values of the mesh $u(X_{d})$ are periodic in
time and $u(X_{d})$ is a discrete travelling wave, i.e.
$u^{i}_{j}=u^{i+zk_{0}}_{j+zk_{1}}$ for all $z\in\mathbb{Z}$, where
$\frac{c\Delta t}{\Delta x}=\frac{k_{0}}{k_{1}}$ with $k_{0}\in\mathbb{Z}$,
$k_{1}\in\mathbb{N}$. The theory of lattice differential equationsMallet-Paret
(2003); Hupkes and Vleck (2022) provides a way to treat $u(X_{d})$ as a
travelling wave even if $\frac{c\Delta t}{\Delta x}\not\in\mathbb{Q}$.

To design machine learning architectures for discrete field theories that
effectively capture solutions corresponding to travelling waves in the
continuous theory, we pay close attention to invariances of the field theory
that relate to structural properties of the dynamical system that governs
travelling waves. In particular, autonomy of a field theory relates to the
fact that travelling waves themselves are governed by a variational principle.
This is made precise in the following Section.

######  Proposition 3 (Symmetric criticality for travelling waves).

Let $c\in\mathbb{R}\setminus\\{0\\}$, $b>0$,
$X=([0,b/c]/\sim)\times([0,b]/\sim)$ be a torus (periodic boundary
conditions), $W=C^{1}(X,\mathbb{R}^{d})$ continuously differentiable functions
and $\Sigma=\\{u\in W\,|\,u(t+s,x+cs)=u(t,x)\\}$ be travelling waves with wave
speed $c$. Let $S\colon W\to\mathbb{R}$

$S(u)=\int_{0}^{b/c}\int_{0}^{b}L(u(t,x),u_{t}(t,x),u_{x}(t,x))\mathrm{d}x\mathrm{d}t$
(40)

be a continuously differentiable action functional with an autonomous
Lagrangian $L\colon(\mathbb{R}^{d})^{3}\to\mathbb{R}$. Denote the restriction
of $S$ to $\Sigma$ by $S_{\Sigma}\colon\Sigma\to\mathbb{R}$. Then a travelling
wave $u\in\Sigma$ is a stationary point of $S$ if and only if $u$ is a
stationary point of $S_{\Sigma}$. Using the identification $\Sigma\cong
C^{1}([0,b]/\sim,\mathbb{R}^{d})$, $u(t,x)=f(x-ct)=f(\xi)$ for $u\in\Sigma$,
we have

$S_{\Sigma}(f)=\frac{b}{c}\int_{0}^{b}L(f(\xi),-cf_{\xi}(\xi),f_{\xi}(\xi))\mathrm{d}\xi.$
(41)

######  Remark 14.

The statement of 3 holds true for standing travelling waves ($c=0$), when the
fraction $b/c$ in the expression of the domain $X$ as well as in (40) and in
(41) are substituted by 1. To include the case $c=0$ we will omit the
irrelevant factor $b/c$ when reporting $S_{\Sigma}$.

Notice that for $u\in\Sigma$ to be a stationary point of $S_{\Sigma}$, the
action needs to be stationary at $u$ with respect to variations through
travelling waves. To be a stationary point of $S$, however, $u$ needs to be
stationary with respect to all variations $\delta u\colon X\to\mathbb{R}^{d}$.
The Section says that stationarity of $S_{\Sigma}$ is already sufficient to
conclude stationarity of $S$.

The proof checks invariance of $S(u)$ under the transformation $u\mapsto s.u$
with $s.u(t,x)=u(t+s,x+cs)$ and verifies the conditions of Palais’ principle
of symmetric criticality (PSC)Palais (1979); Torre (2011). A proof is
contained in Appendix C.

In particular, 3 says that travelling waves are again governed by a
variational principle, namely by the restricted functional
$S_{\Sigma}\colon\Sigma\to\mathbb{R}$. We illustrate this on the wave equation
and the Schrödinger equation.

######  Example 9 (PSC for wave equation).

Consider the wave equation (2) with periodic boundary conditions. Let
$\xi=x-ct$. By 3, travelling wave solutions $u(t,x)=f(x-ct)$ of (4) exactly
correspond to the stationary points of
$S_{\Sigma}(f)\colon[0,b]/\sim\to\mathbb{R}$ with

$\begin{split}S_{\Sigma}=\int_{0}^{b}\left(\frac{1}{2}(c^{2}-1)f_{\xi}^{2}-V(f)\right)\mathrm{d}\xi.\end{split}$
(42)

Although 3 and 14 guarantee variational structure of travelling wave solutions
for any $c\in\mathbb{R}$, they do not imply existence of non-trivial
solutions: for $V(f)=\frac{1}{2}f^{2}$ a Fourier ansatz in the Euler–Lagrange
equation $(1-c^{2})f_{\xi\xi}=\nabla V(f)$ reveals the condition

$c^{2}=1+\frac{l^{2}}{4\pi^{2}m^{2}},\quad m\in\mathbb{Z}\setminus\\{0\\}$

with solutions
$u(t,x)=\alpha_{1}\sin(\kappa(x-ct))+\alpha_{2}\cos(\kappa(x-ct))$, where
$\kappa=2\pi m/b$, $\alpha_{1},\alpha_{2}\in\mathbb{R}$.

######  Example 10 (PSC in Schrödinger equation).

In analogy to 9, travelling waves $\Psi(t,x)=f(x-ct)=f(\xi)\in\mathbb{C}$ with
wave speed $c\in\mathbb{R}$ in the Schrödinger equation (4) with periodic
boundary conditions are governed by the variational principle
$S_{\Sigma}\colon[0,b]/\sim\to\mathbb{R}$ with
$S_{\Sigma}(f)=\int_{0}^{b}L_{\Sigma}(f,f_{\xi})\mathrm{d}\xi$ with

$L_{\Sigma}(f,f_{\xi})=\frac{c\hbar}{2}(\overline{f}f_{\xi}-f\overline{f_{\xi}})-f_{\xi}\overline{f_{\xi}}-V(f\overline{f}).$
(43)

A Fourier ansatz for $f$ reveals that for any
$m\in\mathbb{Z}\setminus\\{0\\}$, $\alpha\in\mathbb{R}$ there is a travelling
wave $\Psi(t,x)=f(\xi)=\alpha e^{\mathrm{i}\frac{2\pi m}{b}\xi}$ with
$\xi=x-ct$ and

$c=\frac{2\pi m}{\hbar b}-\frac{b}{2\pi\hbar}\nabla V(\alpha^{2}).$ (44)

Let us now obtain corresponding statements for discrete field theories. To
make sense of travelling waves in discrete equations such as (14) and (15),
standard approaches in lattice differential equationsMallet-Paret (2003);
Hupkes and Vleck (2022) replace discrete equations by functional equations. We
recover this interpretation from a variational perspective.

######  Lemma 1 (Continuous action for discrete field theory).

Let $X=([0,b_{0}]/\sim)\times([0,b_{1}]/\sim)$ be a torus and $L_{d}$ the
3-point discrete Lagrangian density of 6. Let $\Delta t,\Delta x>0$. Consider
the action $S\colon X\to\mathbb{R}$,

$S(u)=\int_{X}L_{d}(u(t,x),u(t+\Delta t,x),u(t,x+\Delta
x))\mathrm{d}x\mathrm{d}t.$ (45)

Stationary points $u\colon X\to\mathbb{R}^{d}$ are the solutions of the
functional equation

$\begin{split}\frac{\partial}{\partial
u}\big{(}&\phantom{+}L_{d}(u(t,x),u(t+\Delta t,x),u(t,x+\Delta x))\\\
&+L_{d}(u(t-\Delta t,x),u(t,x),u(t-\Delta t,x+\Delta x))\\\
&+L_{d}(u(t,x-\Delta x),u(t+\Delta t,x-\Delta x),u(t,x))\big{)}=0.\end{split}$
(46)

Notice that (46) is (14) with $u^{i+s}_{j+r}$ substituted by $u(t+s\Delta
t,x+r\Delta x)$. The proof follows by a standard variational calculus.
Analogous statements hold for all autonomous discrete Lagrangians (11), in
particular for the 4-point discrete Lagrangian of 6. In analogy to 3 we have
the following Section.

######  Proposition 4 (PSC for TW in discrete theory).

Let $c\in\mathbb{R}\setminus\\{0\\}$, $b>0$,
$X=([0,b/c]/\sim)\times([0,b]/\sim)$ be a torus, $W=C^{1}(X,\mathbb{R}^{d})$
continuously differentiable functions and $\Sigma=\\{u\in
W\,|\,u(t+s,x+cs)=u(t,x)\\}$ be travelling waves with wave speed $c$. Let
$S\colon W\to\mathbb{R}$ be defined as in (45). Denote the restriction of $S$
to $\Sigma$ by $S_{\Sigma}\colon\Sigma\to\mathbb{R}$. Then a travelling wave
$u\in\Sigma$ is a stationary point of $S$, if and only if $u$ is a stationary
point of $S_{\Sigma}$.

The proof is analogous to the proof of 3 since the computation in (60) of
Appendix C works for all autonomous Lagrangians. An analogous statement holds
in the case of standing travelling waves $c=0$.

######  Corollary 1 (TW in discrete wave equation).

Travelling waves of the discrete wave equation (7) with periodic boundary
conditions in space (period $b$) are governed by
$S_{\Sigma}\colon\Sigma\to\mathbb{R}$ with

$\begin{split}S_{\Sigma}(f)&=\int_{0}^{b}\Bigg{(}\frac{1}{2}\left(\frac{f(\xi+\Delta
t)-f(\xi)}{\Delta t}\right)^{2}\\\
&\phantom{\frac{b}{c}\int_{0}^{b}}-\frac{1}{2}\left(\frac{f(\xi+\Delta
x)-f(\xi)}{\Delta x}\right)^{2}-V(f(\xi))\Bigg{)}\mathrm{d}\xi\end{split}$

using again the identification $\Sigma\cong C^{1}([0,b]/\sim,\mathbb{R})$.

######  Corollary 2 (TW in discrete Schrödinger).

Travelling waves of the discrete Schrödinger equation (8) are governed by
$S_{\Sigma}\colon\Sigma\to\mathbb{R}$ with

$S_{\Sigma}(f)=\int_{0}^{b}L_{d}(f(\xi-\Delta t-\Delta x),f(\xi-\Delta
t),f(\xi),f(\xi-\Delta x))\mathrm{d}\xi$

and with $L_{d}$ of (19), using again the identification $\Sigma\cong
C^{1}([0,b]/\sim,\mathbb{R})$.

The expressions $S_{\Sigma}(f)$ follows from a calculation analogous to (60).
In the case $c\not=0$ we have omitted the irrelevant constant factor
$\frac{b}{c}$ when reporting $S_{\Sigma}$.

######  Remark 15 (Blended backward error analysis).

Using the technique of blended backward error analysisMcLachlan and Offen
(2023, 2022) we show that the discrete action $S_{\Sigma}$ can formally be
written as a continuous first order theory
$S_{\Sigma}(f)=\int\tilde{L}(f,f_{\xi})\mathrm{d}\xi$, where $\tilde{L}$ is a
formal power series in $\Delta t$, $\Delta x$.

Figure 6: First two travelling waves ($m=1,2$ in 16) of the discrete wave
equation with quadratic potential. The wave crest appears rough due to the
mesh resolution. This is expected in discrete theories.

######  Remark 16 (TW in discrete wave equation).

Consider $V(u)=\frac{1}{2}u^{2}$ in the discrete wave equation (7) with
periodic boundary conditions in space (period $b$). A Fourier ansatz for $f$
reveals (away from some resonant or degenerate cases) that
$u(t,x)=f(\xi)=\alpha_{1}\sin(\kappa\xi)+\alpha_{2}\cos(\kappa\xi)$,
$\xi=x-ct$, $\kappa=2\pi m/b$, ($m\in\mathbb{Z}\setminus\\{0\\}$,
$\alpha_{1},\alpha_{2}\in\mathbb{R}$) is a solution with wave speed $c$ a real
solution of

$\cos(\kappa c\Delta t)=1-\frac{\Delta t^{2}}{2}+\frac{\Delta t^{2}}{\Delta
x^{2}}(\cos(\kappa\Delta x)-1).$

See Fig. 6 for plots of the first two travelling waves ($m=1,2$).

Figure 7: Left: Travelling waves in 17 of the discrete Schrödinger equation
($m=1$). The wave crest appears rough due to the mesh resolution. Right:
Spurious travelling wave.

######  Remark 17 (TW in discrete Schrödinger).

Consider a linear potential $V(u)=\beta u$ in the discrete Schrödinger
equation (8) with periodic boundary conditions in space (period $b$). A
Fourier ansatz for $f$ reveals (away from some resonant or degenerate cases)
that $u(t,x)=f(\xi)=\alpha e^{\mathrm{i}\kappa\xi}$, $\xi=x-ct$, $\kappa=2\pi
m/b$, ($m\in\mathbb{Z}\setminus\\{0\\}$, $\alpha\in\mathbb{R}$) is a solution
with wave speed

$c=\frac{2}{\kappa\Delta t}\left(\arctan\left(\frac{2}{\hbar}\frac{\Delta
t}{\Delta x^{2}}\tan^{2}(\frac{1}{2}\kappa\Delta x)+\frac{\beta}{2\hbar}\Delta
t\right)+s\pi\right),$

where $s\in\mathbb{Z}$. Asymptotic expansion in $\Delta t,\Delta x$ recovers
the expression (44) up to higher order terms. Additional travelling wave
solutions of the discrete theory are obtained with wave speed
$c=\frac{b(2\tilde{m}+1)}{2m\Delta t}$, $\tilde{m}\in\mathbb{Z}$. These are
regarded as spurious since they do not converge to solutions of the continuous
theory as $\Delta t\to 0$. (See Fig. 7.)

##  VI Numerical Experiments

We report numerical experiments based on learning discrete field theories for
the discrete wave equation (7) and the discrete Schrödinger equation (8) from
data. The discrete field theories with the discrete Lagrangian densities (17),
(19) are regarded as the true theories rather than as approximations of
continuous theories in this context.

The use-case of our machine learning architecture is to learn effective models
of variational dynamical systems whose underlying dynamics is an unknown
variational partial differential equation. The observational data is discrete
and lies on a spacetime lattice. We, therefore, restrict ourselves to explore
models which capture solutions whose discrete Fourier modes exhibit an
exponential decay. This disregards solutions that do not have a limit as
lattice parameters are made asymptotically small. These are sometimes referred
to as spurious solutions in numerical integration theoryHairer, Lubich, and
Wanner (2013) as they do not correspond to solutions of an underlying
continuous theory.

We evaluate the performance of the data-driven models by computing or
extending solutions of the trained models to initial data and checking
consistency against the true models, compare performance to alternative
reduced order modelling approaches (ROMs), and analyse how well travelling
waves are captured. In the experiments, travelling waves will not be part of
the training data but will be discovered using the trained model.

###  VI.1 Wave equation

Consider the space-time domain $X=[0,T]\times[0,b]/\sim$ (periodic boundary
conditions in space), $T=0.5$, $b=1$, discretisation parameters $\Delta
t=0.025$, $\Delta x=0.05$, and quadratic potential $V(u)=\frac{1}{2}u^{2}$.
The discrete lattice is denoted by $X_{d}$ and $M=b/\Delta x=20$ is the number
of spatial interior grid points and $N=T/\Delta t=20$ the number of time
steps.

####  VI.1.1 Training data

Figure 8: Typical elements in the training data set (left) and their absolute
value of amplitudes of Fourier modes averaged over time (right). Bottom:
Absolute values of amplitudes of Fourier modes averaged over time and all
training data.

To generate initial data, $M^{\mathrm{rfft}}=11$ Fourier coefficients
$\gamma_{j}=Me^{-2j^{4}}\eta_{j}$ ($j=0,\ldots,M^{\mathrm{rfft}}$) are
computed where $\eta_{j}$ are samples from a standard normal distribution. An
application of the inverse of the discrete Fourier transform for real input
($\mathtt{irfft}$ of Julia package FFTW.jlFrigo and Johnson (2005)) yields
initial data $(u^{0}_{j})_{j=0}^{M-1}$. Moreover, initial velocities
$(\dot{u}^{0}_{j})_{j=0}^{M-1}$ are drawn from a standard normal distribution.
Then $K=80$ solutions of the field theory are computed (4). Figure 8 displays
two typical solutions and the absolute values of the amplitudes of their
discrete Fourier modes averaged over all time steps. A corresponding average
over all times and all solutions is displayed at the bottom of Fig. 8. Each
solution $U^{(k)}=(u^{i,(k)}_{j})_{0\leq j<M}^{0\leq i\leq N}$ yields
$(N-1)M=380$ tuples
$(u^{i,(k)}_{j},u^{i+1,(k)}_{j},u^{i,(k)}_{j+1},u^{i-1,(k)}_{j},u^{i-1,(k)}_{j+1},u^{i,(k)}_{j-1},u^{i+1,(k)}_{j-1})_{0\leq
j<M}^{0<i<N}$ of data corresponding to the $7$-point stencil. The 7-point
tuples constitute the training data set $\mathcal{D}$. $\mathcal{D}$ contains
30400 tuples of stencil data. For batch trainingMarsland (2014), $\mathcal{D}$
is partitioned into batches of size 10.

####  VI.1.2 Training

The 3 point discrete Lagrangian of 6 is modelled by a neural networks (3
layers, interior layer with 10 nodes, 160 parameters in total, activation
function $\mathtt{tanh}$). Consider $\ell_{\mathrm{data}}$ of (35)

$\ell_{\mathrm{data}}=\sum_{\mathrm{stencil\,data}\in\mathrm{batch\,of\,}\mathcal{D}}|\mathrm{DEL}(L_{d})(\mathrm{stencil\,data})|^{2}$
(47)

with $\mathrm{DEL}(L_{d})$ of (14), and $\ell_{\mathrm{reg}}$ from (36)

$\ell_{\mathrm{reg}}=\alpha\sum\left\|\left(\frac{\partial^{2}L_{d}}{\partial
u^{i}_{j}\partial
u^{i+1}_{j}}({u^{i}_{j}}^{(k)},{u^{i+1}_{j}}^{(k)},{u^{i}_{j+1}}^{(k)})\right)^{-1}\right\|^{2},$
(48)

with normalisation $\alpha=\frac{1}{\mathrm{batch\,size}}=\frac{1}{10}$. The
sum is taken over all stencil data
$(u^{i,(k)}_{j},u^{i+1,(k)}_{j},u^{i,(k)}_{j+1},u^{i-1,(k)}_{j},u^{i-1,(k)}_{j+1},u^{i,(k)}_{j-1},u^{i+1,(k)}_{j-1})$
of a batch.

For each batch of $\mathcal{D}$ one optimisation step with the optimisation
method $\mathtt{adam}$Kingma and Ba (2015) with the default learning
parameters of the Julia/Flux.jlInnes  _et al._ (2018); Innes (2018)
implementation is applied to the parameters of the neural network using the
loss function $\ell=\ell_{\mathrm{data}}+\ell_{\mathrm{reg}}$. This is
repeated 1320 times (number of epochs). For the trained model we have
$\ell_{\mathrm{data}}\approx 8.6\cdot 10^{-8}$ and $\ell_{\mathrm{reg}}\approx
1.4\cdot 10^{-7}$. (Here $\ell_{\mathrm{data}}$, $\ell_{\mathrm{reg}}$ are
evaluated on the full data set $\mathcal{D}$.)

####  VI.1.3 Evaluation

##### Extrapolation

Figure 9: True solution of the model (left) versus computed solution from the
learned model (right). The initial data
$u^{0}_{j}=u^{1}_{j}=\sin(4\pi/Tj\Delta x)$ was not seen during training. Top
line shows predictions until $T=0.5$ (as in training), while the bottom line
shows extrapolation until $20T=10$.

The trained neural network model of a discrete Lagrangian density is used to
compute solutions of the learned field theory initialised from the unseen
initial value $u^{0}_{j}=u^{1}_{j}=\sin(4\pi j\Delta x)$, $j=0,\ldots,M-1$
(Fig. 9). Computations are performed as explained in 4. On the training domain
the absolute difference between a true solution is
$\|U_{\mathrm{true}}-U_{\mathrm{predicted}}\|_{\infty}\approx 0.057$. On the
extrapolated domain with final time $20T$ the error rises to
$\|U_{\mathrm{true}}-U_{\mathrm{predicted}}\|_{\infty}\approx 0.343$ and to
$\approx 1.41$ for $100T$.

##### Travelling waves

The travelling waves of Fig. 6 fulfil the discrete Euler–Lagrange equations
(14) for the trained discrete Lagrangian density $L_{d}$ with a maximal error
of $\max_{i,j}|\mathrm{DEL}(L_{d})^{i}_{j}|<0.004$ ($m=1$, left plot of Fig.
6) or $\max_{i,j}|\mathrm{DEL}(L_{d})^{i}_{j}|<0.033$ ($m=2$, right plot of
Fig. 6), respectively. Therefore, these true travelling waves constitute
solutions of the learned discrete field theory (up to small error). Travelling
waves can be found numerically in the trained model as well: for a guess of a
wave speed $c$ and a periodic function $f\colon[0,b]/\sim\to\mathbb{R}$

$f(\xi)=\sum_{m=-(1+\lfloor\frac{M}{2}\rfloor)}^{1+\lfloor\frac{M}{2}\rfloor}\hat{f}_{|m|}e^{\frac{2\pi\mathrm{i}m}{b}\xi}$
(49)

represented by $1+\lfloor\frac{M}{2}\rfloor$ Fourier coefficients
$\hat{f}_{m}\in\mathbb{C}$, the objective function

$\ell_{\mathrm{wave}}(c,(\hat{f})_{m=0}^{1+\lfloor\frac{M}{2}\rfloor})=\sum_{i=1}^{N-1}\sum_{j=0}^{M-1}\|\mathrm{DEL}(L_{d})^{i}_{j}\|^{2}$

is minimised with $\mathrm{DEL}(L_{d})^{i}_{j}$ from (14) with
$u^{i}_{j}=f(i\Delta t-cj\Delta x)$ for all $i,j$. To avoid finding the
trivial solution, the normalisation condition
$\ell_{\mathrm{wave}}^{\mathrm{unit}}=|\|f\|^{2}_{L^{2}}-1|$ is added to the
objective. Here $\|f\|^{2}_{L^{2}}=\Delta x\sum|\hat{f}_{|m|}|^{2}$ is the
discrete $L_{2}$ norm of $f$.

For the experiment, we initialise the optimisation procedure from the true
Fourier coefficients $(\hat{f}_{m})_{m}$ and wave speed $c$ corresponding to
the true travelling waves of 16 (wave number $m=1,2$) perturbed by random
noise (identically independently normally distributed, standard deviation
$\sigma=0.5$). The located travelling wave $U^{\mathrm{TW}}$ has maximal norm
error $\|U^{\mathrm{TW}}-U^{\mathrm{TW}}_{\mathrm{ref}}\|_{\infty}<0.116$
(wave number $m=1$) or
$\|U^{\mathrm{TW}}-U^{\mathrm{TW}}_{\mathrm{ref}}\|_{\infty}<0.31$ (wave
number $m=2$), respectively. The wave speeds are found up to an error of
$|c^{\mathrm{TW}}-c^{\mathrm{TW}}_{\mathrm{ref}}|<0.0009$ or
$|c^{\mathrm{TW}}-c^{\mathrm{TW}}_{\mathrm{ref}}|<0.0082$, respectively.
Figure 10 compares the contours of the located and the true discrete
travelling wave for the case $m=2$.

Figure 10: Contour plot of true travelling wave (left) with wave number 2 and
the located travelling wave (right).

###  VI.2 Coarse Mesh

Even though training data $(U^{(k)})_{k=1}^{K}$ (see Section VI.1.1) is
available on the fine mesh $X_{d}$ with $M(N+1)$ mesh points, the discrete
model $L_{d}$ can be trained on a coarser mesh. This can be useful to save
computational costs when only low resolution predictions are required.
Consider the mesh $X_{d}^{\mathrm{coarse}}$ for $X$ with $(\Delta
t^{\mathrm{coarse}},\Delta x^{\mathrm{coarse}})=(2\Delta t,2\Delta x)$. We
find $(N-3)*M*K=27200$ tuples
$(u^{i,(k)}_{j},u^{i+2,(k)}_{j},u^{i,(k)}_{j+2},u^{i-2,(k)}_{j},u^{i-2,(k)}_{j+2},u^{i,(k)}_{j-2},u^{i+2,(k)}_{j-2})_{0\leq
j<M}^{0<i<N-1}$ with data over the coarse 7-point stencil in the observational
data. Notice that this makes efficient use also of data over mesh points in
$X_{d}$ that do not lie in $X_{d}^{\mathrm{coarse}}$. Let
$\mathcal{D}^{\mathrm{coarse}}$ be a set of the 27200 training data tuples.
The discrete Lagrangian density $L_{d}^{\mathrm{coarse}}$ is parametrised with
a neural network of the same architecture as before. For training, the network
parameters are randomly initialised. The same expressions for the loss
function (47), (48) and the same batching and optimisation methods are used.
$L_{d}$ is trained for 63800 epochs with final losses
$\ell_{\mathrm{data}}\approx 1.1\cdot 10^{-5}$, $\ell_{\mathrm{reg}}\approx
2.5\cdot 10^{-5}$.

Figure 11: Model for discrete wave equation trained for a coarser grid. Top:
Computation of solution to unseen initial data
$u^{0}_{j}=u^{1}_{j}=\sin(2\pi/Tj\Delta x^{\mathrm{coarse}})$ on extended
domain $20T$ (left: reference, right: data-driven model). Bottom:
Initialisation with exact data $u^{0}_{j}$, $u^{1}_{j}$ of a travelling wave.
(left: reference, right: data-driven model)

Figure 11 shows numerical experiments with the trained coarse model. In the
top row of the figure, the (coarse) 7-point stencil is used to compute a
solution $U_{\mathrm{predicted}}$ on $X_{d}^{\mathrm{coarse}}$ initialised
from the unseen initial data $u^{0}_{j}=u^{1}_{j}=\sin(2\pi/Tj\Delta
x^{\mathrm{coarse}})$. The prediction contains an extrapolation until
$20T=10$. A reference is computed with the 7 point stencil (reducing to a 5
point stencil) w.r.t. the original, fine mesh $X_{d}$ applied to the discrete
Lagrangian density (17) of 7 (the true model). Subsampling over
$X_{d}^{\mathrm{coarse}}$ yields the reference $U_{\mathrm{true}}$. We obtain
$\|U_{\mathrm{predicted}}-U_{\mathrm{true}}\|_{\infty}\approx 0.33$. The
bottom row of Fig. 11 shows the prediction
$U_{\mathrm{predict}}^{\mathrm{TW}}$ of the model when initialised from data
$u^{0}_{j}=u^{1}_{j}$ of a travelling wave $(m=1)$ of the true system. We have
$\|U_{\mathrm{predict}}^{\mathrm{TW}}-U_{\mathrm{true}}^{\mathrm{TW}}\|_{\infty}\approx
0.24$.

###  VI.3 Comparison to MOR based technique

For comparison, we employ model order reduction in the spatial variable and
then learn the temporal dynamics. This approach was explored in a recent
articles in the continuous Lagrangian settingSharma and Kramer (2022) and in
related settings and problemsAllen-Blanchette  _et al._ (2020); Mason  _et
al._ (2022); Buchfink, Glas, and Haasdonk (2023); Carlberg, Tuminaro, and
Boggs (2015); Sharma  _et al._ (2023); Sharma, Wang, and Kramer (2022);
Tyranowski and Kraus (2022). However, for consistency we will continue using
discrete theory and continue using the training data of the discrete wave
equation model (Section VI.1.1).

We use standard principle component analysisBrunton and Kutz (2019); Trefethen
and Bau (1997) to identify a linear projection map
$\mathrm{pr}\colon\mathbb{R}^{M}\to\mathbb{R}^{M_{\mathrm{red}}}=:Q$ and a
recovery map $R\colon Q\to\mathbb{R}^{M}$ such that $\mathrm{pr}\circ
R=\mathrm{id}_{Q}$ with $M_{\mathrm{red}}<M$. Then a neural network model of a
discrete Lagrangian $L^{Q}_{d}\colon Q\times Q\to\mathbb{R}$ is trained to
recover the temporal dynamics. To identify $\mathrm{pr}$, let

$\mathcal{U}^{\mathrm{space}}=\begin{pmatrix}U^{(1)},&U^{(2)},&\ldots,&U^{(K)}\end{pmatrix}\in\mathbb{R}^{MK(N+1)}$

denote a matrix arranging the matrices $U^{(k)}=(u^{i}_{j})_{0\leq i\leq
N}^{0\leq j<M}$ of training data horizontally. Consider its singular value
decomposition $\mathcal{U}^{\mathrm{space}}=A\Sigma B^{\top}$. Let
$A^{\prime}=\begin{pmatrix}a_{1},&a_{2}\end{pmatrix}\in\mathbb{R}^{M\times
M_{\mathrm{red}}}$ contain the first $M_{\mathrm{red}}$ columns of $A$ and set
$\mathrm{pr}(u)={A^{\prime}}^{\top}u$ and $R(q)=A^{\prime}q$. For
$M_{\mathrm{red}}=2$ the average relative reconstruction error is

$\frac{1}{MK(N+1)}\sum_{k=1}^{MK(N+1)}\frac{\|\mathfrak{u}_{k}-R(\mathrm{pr}(\mathfrak{u}_{k}))\|_{2}}{\|\mathfrak{u}_{k}\|_{2}}\approx
0.0026,$

where $\mathfrak{u}_{k}$ denote the columns of $\mathcal{U}^{\mathrm{space}}$.

To learn the temporal dynamics on the latent space, consider a neural network
model of $L^{Q}_{d}\colon Q\times Q\to\mathbb{R}$ (3 layers, interior layer of
size 10, 170 parameters, activation function $\mathtt{softplus}$). Let
$q^{i,(k)}\in Q$ denote the projected training data
$q^{i}_{(k)}=\mathrm{pr}(({u^{i}_{j}}^{(k)})_{0\leq j<M})$ ($1\leq k\leq
K,0\leq i\leq N$) and let

$\begin{split}\mathcal{D}^{Q}=\big{\\{}(q^{i-1}_{(k)},q^{i}_{(k)},q^{i+1}_{(k)})&\,|\,1\leq
k\leq K,0<i<N\big{\\}}\end{split}$

be the collection of tipples in the projected training data. Consider the loss
function $\ell^{Q}=\ell_{\mathrm{data}}^{Q}+\ell_{\mathrm{reg}}^{Q}$ with

$\ell_{\mathrm{data}}^{Q}=\sum_{q_{\mathrm{stencil}}\in\mathcal{D}^{Q}}\mathrm{DEL}(L_{d}^{Q})(q_{\mathrm{stencil}})^{2}$
(50)

with

$\mathrm{DEL}(L_{d}^{Q})(q_{\mathrm{stencil}})=\frac{\partial}{\partial
q^{i}}\left(L_{d}^{Q}(q^{i-1},q^{i})+L_{d}^{Q}(q^{i},q^{i+1})\right)$ (51)

where $q_{\mathrm{stencil}}=(q^{i-1},q^{i},q^{i+1})$ and regularisation

$\ell_{\mathrm{reg}}^{Q}=\sum_{q_{\mathrm{stencil}}\in\mathcal{D}^{Q}}\left\|\left(\frac{\partial^{2}L_{d}^{Q}}{\partial
q^{i}\partial q^{i+1}}L_{d}^{Q}(q^{i},q^{i+1})\right)^{-1}\right\|^{2}.$

The regularisation is an analogy to (36) for discrete Lagrangians of two input
arguments optimised for forward propagation by solving (51) for $q^{i+1}$. The
summands are approximated with 3 inverse vector iterations (see 9). Again,
batch training is performed with a batch size of 10. The regularisation term
$\ell_{\mathrm{reg}}^{Q}$ was scaled by $10^{-8}$ such that both terms of the
loss function are of the same order of magnitude at the beginning of training
with randomly initialised parameters of the neural network. After 54750 epochs
of optimisation with $\mathtt{adam}$, $\ell_{\mathrm{data}}^{Q}<1.19\cdot
10^{-5}$, $\ell_{\mathrm{reg}}^{Q}<1.54\cdot 10^{-5}$.

Figure 12: While the ROM model successfully predicts motions initialised close
to training data (top left) it struggles to generalise (top right compare to
top of Fig. 9). Moreover, initialising with data of a true travelling wave
does not produce a travelling wave (bottom left), while spurious travelling
waves fulfil the discrete Euler–Lagrange equations when projected to the
latent space (bottom right).

The trained model can recover the training data shown in Fig. 8 in max-error
norm $\|\cdot\|_{\infty}$ up to 0.003 or 0.006 when initialised with the exact
data at time $t=0$, $t=\Delta t$. Figure 12 shows that while the ROM model
successfully predicts dynamics with initial data close to training data (top
left), it fails on initial data such as $u^{0}_{j}=u^{1}_{j}=\sin(4\pi j\Delta
x)$, $j=0,\ldots,M-1$ (top right). This can be compared to the proposed
discrete field theory model, which shows good performance on this task (see
Fig. 9).

Unsurprisingly, as the training data does not contain any travelling waves, it
fails to produce a good prediction when initialised with the exact initial
data of a travelling wave. Indeed, the reconstructed solution has no
similarity to a travelling waves (Fig. 12 bottom left). In contrast, the
proposed stencil based scheme produces an output that is visually
undistinguishable (not displayed) from the exact travelling wave (left of Fig.
6).

The travelling wave locator for discrete field theories described in Section
VI.1.3 can easily be adapted to the ROM setting: as a loss function we use

$\ell_{\mathrm{wave}}^{Q}(c,(\hat{f})_{m=0}^{1+\lfloor\frac{M}{2}\rfloor})=\sum_{i=1}^{N-1}\left\|\frac{\partial}{\partial
q^{i}}\left(L_{d}^{Q}(q^{i-1},q^{i})+L_{d}^{Q}(q^{i},q^{i+1})\right)\right\|^{2}$

with $q^{i}=\mathrm{pr}((u^{i}_{j})_{0\leq j<M})$, with $u^{i}_{j}=f(i\Delta
t-cj\Delta x)$ for all $i,j$, and $f$ as in (49). We obtain the result
$U^{\mathrm{TW}}_{Q}$ to the bottom right of Fig. 12: indeed, the projection
$\mathrm{pr}(U^{\mathrm{TW}}_{Q})$ of $U^{\mathrm{TW}}_{Q}$ fulfils the
discrete Euler-Lagrange equation for $L_{d}^{Q}$ on the latent space $Q$ up to

$\left|\sum_{q_{\mathrm{stencil}}\in\mathcal{D}^{Q}_{\mathrm{TW}}}\mathrm{DEL}(L_{d}^{Q})(q_{\mathrm{stencil}})^{2}\
\right|<3\cdot 10^{-9},$ (52)

where $\mathcal{D}^{Q}_{\mathrm{TW}}$ denotes the triples in
$\mathrm{pr}(U^{\mathrm{TW}}_{Q})$. However, it is apparently not an accurate
approximation of the exact travelling wave (left of Fig. 6).

We conclude that the ROM-based method achieves accurate predictions for
initial values close to its training data but can struggles in generalisation
tasks as well as in the travelling wave task, in contrast to the proposed
stencil based model. This is expected since the latent space was identified
based on training data which does not contain travelling waves.

###  VI.4 Schrödinger equation

Consider the space-time domain $X=[0,T]\times[0,b]/\sim$, $T=0.12$, $b=1$ with
periodic boundary conditions in space, discretisation parameters $\Delta
t=0.01$, $\Delta x=0.125$, and linear potential $V(r)=r$. The discrete lattice
is denoted by $X_{d}$ and $M=1/\Delta x=8$ is the number of spatial interior
grid points and $N=T/\Delta t=12$ the number of time steps.

To obtain training data, we use the 4-point Lagrangian density of 8 and
compute 80 solutions: identifying $\mathbb{C}\cong\mathbb{R}^{2}$,
$\Psi=\phi+\mathrm{i}p\mapsto u=\begin{pmatrix}\phi,&p\end{pmatrix}^{\top}$,
the corresponding Lagrangian $L_{\Delta x}^{\Delta t}$ of (25) is degenerate
and of the form as in 3. Therefore, position data at time $t=0$ is sufficient
to initialise the numerical scheme of 4 and no velocity data is required.
Initial position data is generated based on randomly sampled, exponentially
decaying Fourier modes in analogy to Section VI.1.1. Figure 13 (top left)
shows the real part of such a solution.

Figure 13: Learned discrete Schrödinger equation. Top left: Sample of solution
in training data (real part). Top right: Learned model reproduces travelling
wave when initialised from exact travelling wave data (compare to left of Fig.
7). Bottom: Extrapolation of the solution with true discrete Lagrangian (left)
and with the learned discrete Lagrangian (right) (error 0.01 in infinity
norm).

A 4-point Lagrangian $L_{d}\colon(\mathbb{R}^{2})^{4}\to\mathbb{R}$ (see 6) is
modelled as a neural network (2 hidden layers of width 12, activation
$\mathtt{softplus}$) and fitted to the training data. As a loss function
$\ell=\ell_{\mathrm{data}}+\ell_{\mathrm{reg}}$ we adapt the data consistency
term $\ell_{\mathrm{data}}$ of (47) to 4-point Lagrangians. As regularisation
term $\ell_{\mathrm{reg}}$ we use the tamed regulariser (39) to obtain a field
theory optimised for temporal forward propagation. A loss of
$\ell_{\mathrm{data}}<8.04\cdot 10^{-5}$ and $\ell_{\mathrm{reg}}=0$ is
reached after 1544 epochs of batch learning, where a batch consists of 2
blocks and a block contains all stencil data of 3 consecutive time steps.

Figure 13 shows that the trained model successfully extrapolates training data
with an extrapolation error
$\||\Psi_{\mathrm{true}}-\Psi_{\mathrm{predicted}}|_{\mathbb{C}}\|_{\infty}\approx
0.01$ with $T^{\mathrm{extrapolate}}=50\Delta t>12\Delta t=T$ (bottom
figures). Here we have used the notation
$\||\Phi|_{\mathbb{C}}\|_{\infty}=\max_{i,j}|\Phi^{i}_{j}|_{\mathbb{C}}$ for a
complex matrix $\Phi\in\mathbb{C}^{(N+1)\times M}$.

Moreover, when the learned model is initialised with data
$u^{0}=(u^{0}_{j})_{0\leq j<M}$, $u^{1}=(u^{1}_{j})_{0\leq j<M}$ of an exact
travelling wave (17, $m=1$), then it reproduces the exact travelling wave up
to an error of
$\||\Psi_{\mathrm{true}}-\Psi_{\mathrm{predicted}}|_{\mathbb{C}}\|_{\infty}\approx
0.19$. (Compare top right of Fig. 13 showing $\Psi_{\mathrm{predict}}$ with
the left of Fig. 7 showing $\Psi_{\mathrm{true}}$.) Notice that travelling
wave solutions are captured well in the trained model although the training
data does not contain any travelling waves.

##  VII Summary and Future Work

We have developed a machine learning framework to learn discrete field
theories from discrete observation data. Neural network models of discrete
Lagrangian densities are fitted such that the discrete Euler Lagrange
equations are fulfilled on the training data. The discrete Lagrangian
densities are based on local stencils motivated by the theory of variational
integration methodsMarsden and West (2001). Thus, the local nature of
differential equations is reflected in the structure of our approach.

As Lagrangian densities are not uniquely determined by a dynamical system’s
motions, we present a way to systematically derive regularisation terms for
the training procedure. The regularisers guide the optimisation process to a
numerical analysis informed model of a discrete Lagrangian density with good
convergence properties when used in numerical simulations.

In comparison to architectures that project spatial dimensions to a data-
driven latent space (data-driven reduced order models - ROMs) and then learn
the temporal dynamics, our framework learns the discrete field theory on a
space-time lattice. A reduction of computational complexity of the learning
task in our approach is obtained by the use of local stencils rather than by a
contraction of dimensionality. By a comparison with a simple, data-driven ROM
based method, we demonstrate that learning on the lattice can be advantageous
when highly symmetric solutions such as travelling waves are sought in the
data-driven model but are not present in the training data set. Moreover,
locality, symmetries, and fundamental principles such as Palais’ principle of
symmetric criticality can be easier to preserve.

However, while locality of the stencil can be exploited during training in the
evaluation of a data consistency term $\ell_{\mathrm{data}}$, depending on
boundary conditions, numerical stability concerns, or particularities of an
employed stencil, numerical integration of the trained model of a discrete
field theory can require solving high-dimensional problems. Moreover, when the
regulariser employed during training is based on the mentioned numerical
method, then it can inherit the computational complexity. While we have shown
that one can train the field theory on a coarse meshes while efficiently using
the observed training data from the finer mesh, another option is to develop
machine learning architectures tailored to numerical integrators developed for
high-dimensional problems such as low-rank integratorsWalach (2019); Lubich,
Vandereycken, and Walach (2018). This will be the subject of future work.

Further, the machine learning architecture based on discrete Lagrangian
densities for local stencils learns structurally simple solutions well and
extrapolates correctly to highly symmetric solutions, such as travelling
waves, even when these are not in the training data set. Numerical integration
theory offers a diverse set of discretisation methods suitable for a variety
of purposes including handling of shock waves or highly oscillatory solutions.
In future work, we plan to exploit these methods within the machine learning
context to develop reliable architectures for systems with challenging
dynamical features.

Additionally, we propose to learn a representation of a field theory that is
as symmetric as possible to improve qualitative features of numerical
simulations as we have done in the ode caseLishkova  _et al._ (2023). A
combination with this technique will allow for the identification of
conservation laws of the underlying dynamical system.

## Data availability

The source code of the numerical experiments can be found at
https://github.com/Christian-Offen/DLNN_pde.

###### Acknowledgements.

C. Offen acknowledges the Ministerium für Kultur und Wissenschaft des Landes
Nordrhein-Westfalen and computing time provided by the Paderborn Center for
Parallel Computing (PC2).

##  Appendix A Exact discrete Lagrangians

We can view the various discrete Lagrangians introduced in Section II.2 as
discretisations of exact discrete Lagrangians. An exact discrete Lagrangian to
a continuous field theory has the property that solutions to the discrete
Euler–Lagrange equations coincide with the motions of the continuous theory on
the grid. For this we extend our definition of discrete Lagrangians from
Section II.2. For a more detailed introduction and further information we
refer to Marsden  _et al._ (2001).

With the notation of Section II.2, let $\partial X^{l}$ denote the boundary of
the hypercube $X^{l}$. Let $\partial X_{d}=\bigcup_{l\in\mathcal{I}}\partial
X^{l}\subset\mathbb{R}^{n+1}$ denote the $n$-skeleton of $X$.

######  Definition 1 (Discrete Lagrangian and action).

A discrete Lagrangian $L_{d}\colon\bigcup_{l}(\partial
X^{l}\times\mathbb{R}^{k})\to\mathbb{R}$ maps any graph of a boundary function
$u_{\partial X^{l}}\colon\partial X^{l}\to\mathbb{R}^{k}$ to a real number.
The corresponding discrete action $S_{d}\colon\partial
X_{d}\times\mathbb{R}^{k}\to\mathbb{R}$ maps any graph of a boundary function
$u_{\partial X_{d}}\colon\partial X_{d}\to\mathbb{R}^{k}$ on the $n$-skeleton
to

$S_{d}(u_{\partial X_{d}})=\sum_{l\in\mathcal{I}}L_{d}(u_{\partial
X_{d}}|_{\partial X^{l}}).$ (53)

Here $u_{\partial X_{d}}|_{\partial X^{l}}$ denotes the restriction of
$u_{\partial X_{d}}$ to the boundary $\partial X^{l}$ of the cube $X^{l}$.

Notice that in 1 $L_{d}$ is allowed to depend on the graph of the boundary
function rather than just its values, which is implicit in the notation.

######  Definition 2 (Discrete Euler–Lagrange equations).

Let $L_{d}$ be a discrete Lagrangian. For each interior vertex
$v\in\mathring{X}_{d}$ let $\mathcal{I}_{v}\subset\mathcal{I}$ contain the
indices of hypercubes $X^{l}$ containing $v$. The discrete Euler–Lagrange
equations are given as

$\sum_{l\in I_{v}}\frac{\delta L_{d}}{\delta u_{\partial X^{l}}}(u_{\partial
X_{d}}|_{\partial X^{l}})=0$ (54)

for $u_{\partial X_{d}}\colon\partial X_{d}\to\mathbb{R}^{k}$. Here
$\frac{\delta L_{d}}{\delta u_{\partial X^{l}}}$ denotes the variational
derivative.

######  Example 11 (Exact discrete Lagrangian).

For a Lagrangian $L$ as in (1) we define the exact discrete Lagrangian

$L_{d}^{\mathrm{exact}}(u_{\partial
X^{l}})=\int_{X^{l}}L(x,u(x),u_{x_{0}}(x),\ldots,u_{x_{n}}(x))\mathrm{d}x,$
(55)

where $u\colon X^{l}\to\mathbb{R}^{k}$ is the solution (assuming
existence/uniqueness) of the Euler-Lagrange equation $\mathrm{EL}(L)=0$
subject to the boundary condition $u|_{\partial X^{l}}=u_{\partial X^{l}}$.

We conclude with the following observations.

  * •

If $u\colon X\to\mathbb{R}^{k}$ solves the Euler–Lagrange equation (2) then
$S(u)=S_{d}^{\mathrm{exact}}(u|_{\partial X_{d}})$, when $S_{d}$ is the
discrete action (53) for the exact discrete Lagrangian
$L_{d}=L_{d}^{\mathrm{exact}}$.

  * •

In the ode case $\dim X=1$ we have $\partial
X_{d}=X_{d}=\\{x_{0}^{0},\ldots,x_{0}^{N_{0}}\\}$ and the exact discrete
action is defined for functions on the mesh $X_{d}$.

  * •

The discrete Lagrangians (17) of 7 and (19) of 8 can be obtained from (55) by
quadrature.

##  Appendix B Convergence of Newton Iterations to solve the 7-point stencil

We provide a proof of 1.

######  Proposition.

Let $u^{i}_{j}$, $u^{i+1}_{j}$, $u^{i}_{j+1}$, $u^{i-1}_{j}$, $u^{i-1}_{j+1}$,
$u^{i}_{j-1}$, $u^{i+1}_{j-1}$ such that the 7-point stencil (14) is
fulfilled. Let $\mathcal{O}\subset\mathbb{R}^{d}$ be a convex, neighbourhood
of $u^{\ast}=u^{i+1}_{j}$, $\|\cdot\|$ a norm of $\mathbb{R}^{d}$ inducing an
operator norm on $\mathbb{R}^{d\times d}$. Define

$p(u):=\frac{\partial^{2}L_{d}}{\partial u^{i}_{j}\partial
u}(u^{i}_{j},u,u^{i}_{j+1})\in\mathbb{R}^{d\times d}$

and let $\theta$ and $\overline{\theta}$ be Lipschitz constants on
$\mathcal{O}$ for $p$ and for $\mathrm{inv}\circ p$, respectively, where
$\mathrm{inv}$ denotes matrix inversion. Let

$\rho^{\ast}:=\left\|\mathrm{inv}(p(u^{\ast}))\right\|=\left\|\left(\frac{\partial^{2}L_{d}}{\partial
u^{i}_{j}\partial
u^{\ast}}(u^{i}_{j},u^{\ast},u^{i}_{j+1})\right)^{-1}\right\|$ (56)

and let $f(u^{(n)})$ denote the left hand side of (14) with $u_{j}^{i+1}$
replaced by $u^{(n)}$. If
$\|{u}^{(0)}-{u}^{\ast}\|\leq\min\left(\frac{\rho^{\ast}}{\overline{\theta}},\frac{1}{2\theta\rho^{\ast}}\right)$
for ${u}^{(0)}\in\mathcal{O}$, then the Newton Iterations

${u}^{(n+1)}:={u}^{(n)}-\mathrm{inv}(p(u^{(n)}))f(u^{(n)})$

converge quadratically against ${u}^{\ast}$, i.e.

$\|{u}^{(n+1)}-{u}^{\ast}\|\leq\rho^{\ast}\theta\|{u}^{(n)}-{u}^{\ast}\|^{2}.$
(57)

###### Proof.

We adapt standard estimates for Newton’s method (see (Deuflhard and Hohmann,
2003b, §4), for instance) to the considered setting. Let
$f\colon\mathcal{O}\to\mathbb{R}^{d}$ with

$\begin{split}f(u)=\frac{\partial}{\partial
u^{i}_{j}}\big{(}&\phantom{+}L_{d}(u^{i}_{j},u,u^{i}_{j+1})+L_{d}(u^{i-1}_{j},u^{i}_{j},u^{i-1}_{j+1})\\\
&+L_{d}(u^{i}_{j-1},u^{i+1}_{j-1},u^{i}_{j})\big{)}\end{split}$

With this definition, $f(u^{\ast})=0$, $\theta$ is a Lipschitz constant for
$\mathrm{D}f\colon\mathcal{O}\to\mathbb{R}^{d\times d}$,
$u\mapsto\mathrm{D}f(u)$, $\overline{\theta}$ is a Lipschitz constant for
$\mathrm{inv}\circ\mathrm{D}f\colon\mathcal{O}\to\mathbb{R}^{d\times d}$,
$u\mapsto\mathrm{D}f(u)^{-1}$, and
$\rho^{\ast}=\|\mathrm{D}f(u^{\ast})^{-1}\|$. Here $\mathrm{D}f(u)$ denotes
the Jacobian matrix of $f$ at $u\in\mathcal{O}$.

Assume that for $n\in\mathbb{N}$ an iterate $u^{(n)}\in\mathcal{O}$ fulfils
$\|u^{(n)}-u^{\ast}\|\leq\min\left(\frac{\rho^{\ast}}{\overline{\theta}},\frac{1}{2\theta\rho^{\ast}}\right)$.
Then

$\begin{split}\|\mathrm{D}f(u^{(n)})^{-1}&\|=\|\mathrm{D}f(u^{(n)})^{-1}-\mathrm{D}f(u^{\ast})^{-1}+\mathrm{D}f(u^{\ast})^{-1}\|\\\
&\leq\|\mathrm{D}f(u^{(n)})^{-1}-\mathrm{D}f(u^{\ast})^{-1}\|+\rho^{\ast}\\\
&\leq\underbrace{\overline{\theta}\|u^{(n)}-u^{\ast}\|}_{\leq\rho^{\ast}}+\rho^{\ast}\leq
2\rho^{\ast}.\end{split}$

For the next iterate $u^{(n+1)}=u^{(n)}-\mathrm{D}f(u^{(n)})^{-1}f(u^{(n)})$
the distance to $u^{\ast}$ can be bounded:

$\displaystyle\|$ $\displaystyle u^{(n+1)}-u^{\ast}\|$
$\displaystyle=\|u^{(n)}-u^{\ast}-\mathrm{D}f(u^{(n)})^{-1}(f(u^{(n)})-\underbrace{f(u^{\ast})}_{=0})\|$
$\displaystyle\leq\underbrace{\|\mathrm{D}f(u^{(n)})^{-1}\|}_{\leq
2\rho^{\ast}}\Big{\|}\underbrace{\mathrm{D}f(u^{(n)})(u^{(n)}-u^{\ast})}_{=\int_{0}^{1}\mathrm{D}f(u^{(n)})(u^{(n)}-u^{\ast})\mathrm{d}t}$
$\displaystyle\phantom{\|\mathrm{D}f(u^{(n)})^{-1}\|}-\underbrace{(f(u^{(n)})-f(u^{\ast}))}_{\int_{0}^{1}\mathrm{D}f(u^{(n)}+t(u^{\ast}-u^{(n)}))(u^{(n)}-u^{\ast})\mathrm{d}t}\Big{\|}$
$\displaystyle\leq
2\rho^{\ast}\Big{\|}\int_{0}^{1}\left(\mathrm{D}f(u^{(n)})-\mathrm{D}f(u^{(n)}+t(u^{\ast}-u^{(n)}))\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot(u^{(n)}-u^{\ast})\mathrm{d}t\Big{\|}$
$\displaystyle\leq
2\rho^{\ast}\int_{0}^{1}\underbrace{\left\|\mathrm{D}f(u^{(n)})-\mathrm{D}f(u^{(n)}+t(u^{\ast}-u^{(n)}))\right\|}_{\leq\theta
t\|u^{(n)}-u^{\ast}\|}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\left\|u^{(n)}-u^{\ast}\right\|\mathrm{d}t$
$\displaystyle\leq
2\rho^{\ast}\theta\left\|u^{(n)}-u^{\ast}\right\|^{2}\int_{0}^{1}t\mathrm{d}t$
$\displaystyle=\rho^{\ast}\theta\left\|u^{(n)}-u^{\ast}\right\|^{2}$

Moreover, since
$\left\|u^{(n)}-u^{\ast}\right\|\leq(2\theta\rho^{\ast})^{-1}$, we have
$\|u^{(n+1)}-u^{\ast}\|\leq\frac{1}{2}\|u^{(n)}-u^{\ast}\|$. By induction, the
Newton Iterations converge against $u^{\ast}$ and (57) holds true. ∎

##  Appendix C Symmetric criticality of travelling waves in the wave equation

We provide a proof of 3.

######  Proposition (PSC for travelling waves).

Let $c\in\mathbb{R}\setminus\\{0\\}$, $b>0$, $X=[0,b/c]/\sim\times[0,b]/\sim$
be a torus, $W=C^{1}(X,\mathbb{R}^{d})$ continuously differentiable functions
and $\Sigma=\\{u\in W\,|\,u(t+s,x+cs)=u(t,x)\\}$ be travelling waves with wave
speed $c$. Let $S\colon W\to\mathbb{R}$

$S(u)=\int_{0}^{b/c}\int_{0}^{b}L(u(t,x),u_{t}(t,x),u_{x}(t,x))\mathrm{d}x\mathrm{d}t$
(58)

be a continuously differentiable action functional with an autonomous
Lagrangian $L\colon(\mathbb{R}^{d})^{3}\to\mathbb{R}$. Denote the restriction
of $S$ to $\Sigma$ by $S_{\Sigma}\colon\Sigma\to\mathbb{R}$. Then a travelling
wave $u\in\Sigma$ is a stationary point of $S$, if and only if $u$ is a
stationary point of $S_{\Sigma}$. Using the identification $\Sigma\cong
C^{1}([0,b]/\sim,\mathbb{R}^{d})$, $u(t,x)=f(x-ct)=f(\xi)$ for $u\in\Sigma$,
we have

$S_{\Sigma}(f)=\frac{b}{c}\int_{0}^{b}L(f(\xi),-cf_{\xi}(\xi),f_{\xi}(\xi))\mathrm{d}\xi.$
(59)

###### Proof.

Equipped with the norm $\|\cdot\|_{W}$ with $\|u\|_{W}=\max_{(t,x)\in
X}|u(t,x)|+\max_{(t,x)\in X}\|Du(t,x)\|_{\mathbb{R}^{d}\times\mathbb{R}^{2}}$
the space $W$ is a Banach space since the torus $X$ is compact. The group
$([0,b]/\sim,+)$ is compact and acts on the Banach space $W$ by $u\mapsto s.u$
with $s.u(t,x)=u(t+s,x+cs)$. The fixpoints of the group action are $\Sigma$.
The considered action $S$ is invariant under the action, i.e. $S(u)=S(s.u)$
for all $s\in[0,b]/\sim$ since the Lagrangian is autonomous: let $u\in W$,
$s\in[0,b)$.

$\begin{split}S(s.u)&=\int_{0}^{b/c}\int_{0}^{b}L(s.u(t,x),(s.u)_{t}(t,x),(s.u)_{x}(t,x))\mathrm{d}x\mathrm{d}t\\\
&=\int_{-s/c}^{(b-s)/c}\int_{-s}^{b-s}L(u(t,x),u_{t}(t,x),u_{x}(t,x))\mathrm{d}x\mathrm{d}t\\\
&=\int_{0}^{b/c}\int_{0}^{b}L(u(t,x),u_{t}(t,x),u_{x}(t,x))\mathrm{d}x\mathrm{d}t\\\
&=S(u).\end{split}$ (60)

By compactness of the group $([0,b]/\sim,+)$, Palais’ principle of symmetric
criticalityPalais (1979) guarantees that critical points of the restricted
functional $S_{\Sigma}\colon\Sigma\to\mathbb{R}$ constitute critical points of
$S$.

To proof validity of (59), consider coordinates $(x-ct,x+ct)=(\xi,\theta)$ and
a travelling wave $u(t,x)=f(x-ct)=f(\xi)$,
$f\colon[0,b]/\sim\to\mathbb{R}^{d}$. We have

$\begin{split}S_{\Sigma}(f)&=S(u)=\int_{0}^{b/c}\int_{0}^{b}L(f(\xi),-cf_{\xi}(\xi),f_{\xi}(\xi))\mathrm{d}x\mathrm{d}t\\\
&=\int_{-b}^{b}\int_{|\xi|}^{2b-|\xi|}L(f(\xi),-cf_{\xi}(\xi),f_{\xi}(\xi))\frac{\mathrm{d}\theta\mathrm{d}\xi}{2c}\\\
&\stackrel{{\scriptstyle(\ast)}}{{=}}\int_{0}^{b}\int_{0}^{2b}L(f(\xi),-cf_{\xi}(\xi),f_{\xi}(\xi))\frac{\mathrm{d}\theta\mathrm{d}\xi}{2c}\\\
&=\frac{b}{c}\int_{0}^{b}L(f(\xi),-cf_{\xi}(\xi),f_{\xi}(\xi))\mathrm{d}\xi\end{split},$
(61)

where we have used periodicity $f(\xi+b)=f(\xi)$ in $(\ast)$. ∎

######  Remark 18.

For standing travelling waves $c=0$, the torus $X=[0,1]/\sim\times[0,b]/\sim$
is considered. An analogous statement is obtained with $b/c$ substituted by 1
in (58) and (59).

## References

  * Narendra and Parthasarathy (1990) K. Narendra and K. Parthasarathy, “Identification and control of dynamical systems using neural networks,” IEEE Transactions on Neural Networks 1, 4–27 (1990).
  * Han  _et al._ (2019) J. Han, C. Ma, Z. Ma, and W. E, “Uniformly accurate machine learning-based hydrodynamic models for kinetic equations,” Proceedings of the National Academy of Sciences 116, 21983–21991 (2019).
  * Dierkes  _et al._ (2023) E. Dierkes, C. Offen, S. Ober-Blöbaum, and K. Flaßkamp, “Hamiltonian neural networks with automatic symmetry detection,” Chaos: An Interdisciplinary Journal of Nonlinear Science 33 (2023), 10.1063/5.0142969, 063115.
  * Lishkova  _et al._ (2023) Y. Lishkova, P. Scherer, S. Ridderbusch, M. Jamnik, P. Liò, S. Ober-Blöbaum, and C. Offen, “Discrete Lagrangian neural networks with automatic symmetry discovery (accepted),” in  _Proceedings of IFAC World Congress 2023, Yokohama, Japan (9-14/07/2023)_ (IFAC-PaperOnLine, 2023) 2211.10830 .
  * Wu, Xiao, and Paterson (2018) J.-L. Wu, H. Xiao, and E. Paterson, “Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework,” Phys. Rev. Fluids 3, 074602 (2018).
  * Kadupitiya  _et al._ (2020) J. Kadupitiya, F. Sun, G. Fox, and V. Jadhao, “Machine learning surrogates for molecular dynamics simulations of soft materials,” Journal of Computational Science 42, 101107 (2020).
  * Karniadakis  _et al._ (2021) G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, “Physics-informed machine learning,” Nature Reviews Physics 3, 422–440 (2021).
  * Marsden and Ratiu (1999) J. E. Marsden and T. S. Ratiu,  _Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems_ (Springer New York, New York, NY, 1999).
  * Mansfield (2010) E. L. Mansfield,  _A Practical Guide to the Invariant Calculus_, Cambridge Monographs on Applied and Computational Mathematics (Cambridge University Press, 2010).
  * Cranmer  _et al._ (2020) M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, and S. Ho, “Lagrangian neural networks,”  (2020).
  * Hamzi and Owhadi (2021) B. Hamzi and H. Owhadi, “Learning dynamical systems from data: A simple cross-validation perspective, part i: Parametric kernel flows,” Physica D: Nonlinear Phenomena 421, 132817 (2021).
  * Allen-Blanchette  _et al._ (2020) C. Allen-Blanchette, S. Veer, A. Majumdar, and N. E. Leonard, “LagNetViP: A Lagrangian neural network for video prediction (AAAI 2020 symposium on physics guided ai),”  (2020).
  * Lutter, Ritter, and Peters (2019) M. Lutter, C. Ritter, and J. Peters, “Deep Lagrangian Networks: Using physics as model prior for deep learning,” in  _7th International Conference on Learning Representations (ICLR)_ (ICLR, 2019).
  * Greydanus, Dzamba, and Yosinski (2019) S. Greydanus, M. Dzamba, and J. Yosinski, “Hamiltonian Neural Networks,” in  _Advances in Neural Information Processing Systems_, Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Jin  _et al._ (2020) P. Jin, Z. Zhang, A. Zhu, Y. Tang, and G. E. Karniadakis, “SympNets: Intrinsic structure-preserving symplectic networks for identifying Hamiltonian systems,” Neural Networks 132, 166–179 (2020).
  * Bertalan  _et al._ (2019) T. Bertalan, F. Dietrich, I. Mezić , and I. G. Kevrekidis, “On learning Hamiltonian systems from data,” Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 121107 (2019).
  * Chen  _et al._ (2023) Y. Chen, B. Xu, T. Matsubara, and T. Yaguchi, “Variational principle and variational integrators for neural symplectic forms,” in  _ICML Workshop on New Frontiers in Learning, Control, and Dynamical Systems_ (2023).
  * Qin (2020) H. Qin, “Machine learning and serving of discrete field theories,” Scientific Reports 10 (2020), 10.1038/s41598-020-76301-0.
  * Ober-Blöbaum and Offen (2023) S. Ober-Blöbaum and C. Offen, “Variational learning of Euler–Lagrange dynamics from data,” Journal of Computational and Applied Mathematics 421, 114780 (2023).
  * Marsden and West (2001) J. E. Marsden and M. West, “Discrete mechanics and variational integrators,” Acta Numerica 10, 357–514 (2001).
  * Offen and Ober-Blöbaum (2023) C. Offen and S. Ober-Blöbaum, “Learning discrete lagrangians for variational pdes from data and detection of travelling waves,” in  _Geometric Science of Information_, Vol. 14071, edited by F. Nielsen and F. Barbaresco (Springer Nature Switzerland, Cham, 2023) pp. 569–579.
  * Mason  _et al._ (2022) J. Mason, C. Allen-Blanchette, N. Zolman, E. Davison, and N. Leonard, “Learning interpretable dynamics from images of a freely rotating 3d rigid body,”  (2022).
  * Buchfink, Glas, and Haasdonk (2023) P. Buchfink, S. Glas, and B. Haasdonk, “Symplectic model reduction of Hamiltonian systems on nonlinear manifolds and approximation with weakly symplectic autoencoder,” SIAM Journal on Scientific Computing 45, A289–A311 (2023).
  * Carlberg, Tuminaro, and Boggs (2015) K. Carlberg, R. Tuminaro, and P. Boggs, “Preserving Lagrangian structure in nonlinear model reduction with application to structural dynamics,” SIAM Journal on Scientific Computing 37, B153–B184 (2015).
  * Sharma  _et al._ (2023) H. Sharma, H. Mu, P. Buchfink, R. Geelen, S. Glas, and B. Kramer, “Symplectic model reduction of Hamiltonian systems using data-driven quadratic manifolds,” (2023), arXiv:2305.15490 [math.NA] .
  * Sharma and Kramer (2022) H. Sharma and B. Kramer, “Preserving lagrangian structure in data-driven reduced-order modeling of large-scale dynamical systems,” (2022), arXiv:2203.06361 [math.NA] .
  * Sharma, Wang, and Kramer (2022) H. Sharma, Z. Wang, and B. Kramer, “Hamiltonian operator inference: Physics-preserving learning of reduced-order models for canonical Hamiltonian systems,” Physica D: Nonlinear Phenomena 431, 133122 (2022).
  * Tyranowski and Kraus (2022) T. M. Tyranowski and M. Kraus, “Symplectic model reduction methods for the Vlasov equation,” Contributions to Plasma Physics , e202200046 (2022).
  * Palais (1979) R. S. Palais, “The principle of symmetric criticality,” Comm. Math. Phys. 69, 19–30 (1979).
  * Torre (2011) C. G. Torre, “Symmetric Criticality in Classical Field Theory,” AIP Conference Proceedings 1360, 63–74 (2011).
  * Olver (1986) P. J. Olver,  _Applications of Lie Groups to Differential Equations_ (Springer US, 1986).
  * Basdevant (2007) J.-L. Basdevant,  _Variational Principles in Physics_ (Springer New York, New York, NY, 2007).
  * Deriglazov (2009) A. Deriglazov, “On singular Lagrangian underlying the Schrödinger equation,” Physics Letters A 373, 3920–3923 (2009).
  * Marsden  _et al._ (2001) J. E. Marsden, S. Pekarsky, S. Shkoller, and M. West, “Variational methods, multisymplectic geometry and continuum mechanics,” Journal of Geometry and Physics 38, 253–284 (2001).
  * Kieri, Lubich, and Walach (2016) E. Kieri, C. Lubich, and H. Walach, “Discretized dynamical low-rank approximation in the presence of small singular values,” SIAM Journal on Numerical Analysis 54, 1020–1038 (2016).
  * Lubich, Vandereycken, and Walach (2018) C. Lubich, B. Vandereycken, and H. Walach, “Time integration of rank-constrained tucker tensors,” SIAM Journal on Numerical Analysis 56, 1273–1290 (2018), https://doi.org/10.1137/17M1146889 .
  * Walach (2019) H. M. Walach,  _Time integration for the dynamical low-rank approximation of matrices and tensors_ , Ph.D. thesis, Eberhard Karls Universität Tübingen (2019).
  * Schrammer (2022) S. Schrammer,  _On dynamical low-rank integrators for matrix differential equations_ , Ph.D. thesis, Karlsruher Institut für Technologie (KIT) (2022).
  * Rowley and Marsden (2002) C. Rowley and J. Marsden, “Variational integrators for degenerate lagrangians, with application to point vortices,” in  _Proceedings of the 41st IEEE Conference on Decision and Control, 2002._, Vol. 2 (2002) pp. 1521–1527.
  * Kraus (2017) M. Kraus, “Projected variational integrators for degenerate lagrangian systems,” (2017), arXiv:1708.07356 [math.NA] .
  * Ellison  _et al._ (2018) C. L. Ellison, J. M. Finn, J. W. Burby, M. Kraus, H. Qin, and W. M. Tang, “Degenerate variational integrators for magnetic field line flow and guiding center trajectories,” Physics of Plasmas 25 (2018), 10.1063/1.5022277, 052502.
  * Burby, Finn, and Ellison (2022) J. W. Burby, J. M. Finn, and C. L. Ellison, “Improved accuracy in degenerate variational integrators for guiding centre and magnetic field line flow,” Journal of Plasma Physics 88, 835880201 (2022).
  * Deuflhard and Hohmann (2003a) P. Deuflhard and A. Hohmann, “Linear eigenvalue problems,” in  _Numerical Analysis in Modern Scientific Computing: An Introduction_ (Springer New York, New York, NY, 2003) pp. 119–150.
  * Röhrig-Zöllner, Thies, and Basermann (2022) M. Röhrig-Zöllner, J. Thies, and A. Basermann, “Performance of the low-rank TT-SVD for large dense tensors on modern multicore CPUs,” SIAM Journal on Scientific Computing 44, C287–C309 (2022).
  * Marsland (2014) S. Marsland,  _Machine Learning: An Algorithmic Perspective, Second Edition_, 2nd ed. (Chapman & Hall/CRC, 2014).
  * Kingma and Ba (2015) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in  _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_, edited by Y. Bengio and Y. LeCun (2015).
  * Mallet-Paret (2003) J. Mallet-Paret, “Traveling waves in spatially discrete dynamical systems of diffusive type,” in  _Dynamical Systems: Lectures given at the C.I.M.E. Summer School held in Cetraro, Italy, June 19-26, 2000_, edited by J. W. Macki and P. Zecca (Springer Berlin Heidelberg, Berlin, Heidelberg, 2003) pp. 231–298.
  * Hupkes and Vleck (2022) H. J. Hupkes and E. S. V. Vleck, “Travelling waves for adaptive grid discretizations of reaction diffusion systems III: Nonlinear theory,” Journal of Dynamics and Differential Equations (2022), 10.1007/s10884-022-10143-4.
  * McLachlan and Offen (2023) R. I. McLachlan and C. Offen, “Backward error analysis for conjugate symplectic methods,” Journal of Geometric Mechanics 15, 98–115 (2023).
  * McLachlan and Offen (2022) R. I. McLachlan and C. Offen, “Backward error analysis for variational discretisations of pdes,” Journal of Geometric Mechanics 14, 447–471 (2022).
  * Hairer, Lubich, and Wanner (2013) E. Hairer, C. Lubich, and G. Wanner,  _Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations_, Springer Series in Computational Mathematics (Springer Berlin Heidelberg, 2013).
  * Frigo and Johnson (2005) M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proceedings of the IEEE 93, 216–231 (2005), special issue on “Program Generation, Optimization, and Platform Adaptation”.
  * Innes  _et al._ (2018) M. Innes, E. Saba, K. Fischer, D. Gandhi, M. C. Rudilosso, N. M. Joy, T. Karmali, A. Pal, and V. Shah, “Fashionable modelling with Flux,” CoRR abs/1811.01457 (2018), arXiv:1811.01457 .
  * Innes (2018) M. Innes, “Flux: Elegant machine learning with Julia,” Journal of Open Source Software (2018), 10.21105/joss.00602.
  * Brunton and Kutz (2019) S. L. Brunton and J. N. Kutz, “Singular value decomposition (svd),” in  _Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control_ (Cambridge University Press, 2019) p. 3–46.
  * Trefethen and Bau (1997) L. N. Trefethen and D. Bau,  _Numerical Linear Algebra_ (SIAM, 1997).
  * Deuflhard and Hohmann (2003b) P. Deuflhard and A. Hohmann,  _Numerical Analysis in Modern Scientific Computing_ (Springer New York, 2003).