# Discrete Lagrangian Neural Networks with Automatic Symmetry Discovery

Yana Lishkova$\dagger$ Paul Scherer$\dagger$ Steffen Ridderbusch  Mateja
Jamnik  Pietro Liò  Sina Ober-Blöbaum  Christian Offen  University of Oxford,
Oxford, UK  University of Cambridge, Cambridge CB3 0FD, UK  Paderborn
University, D-33098 Paderborn, Germany

###### Abstract

By one of the most fundamental principles in physics, a dynamical system will
exhibit those motions which extremise an action functional. This leads to the
formation of the Euler-Lagrange equations, which serve as a model of how the
system will behave in time. If the dynamics exhibit additional symmetries,
then the motion fulfils additional conservation laws, such as conservation of
energy (time invariance), momentum (translation invariance), or angular
momentum (rotational invariance). To learn a system representation, one could
learn the discrete Euler-Lagrange equations, or alternatively, learn the
discrete Lagrangian function $\mathcal{L}_{d}$ which defines them. Based on
ideas from Lie group theory, in this work we introduce a framework to learn a
discrete Lagrangian along with its symmetry group from discrete observations
of motions and, therefore, identify conserved quantities. The learning process
does not restrict the form of the Lagrangian, does not require velocity or
momentum observations or predictions and incorporates a cost term which
safeguards against unwanted solutions and against potential numerical issues
in forward simulations. The learnt discrete quantities are related to their
continuous analogues using variational backward error analysis and numerical
results demonstrate the improvement such models can have both qualitatively
and quantitatively even in the presence of noise.

###### keywords:

Nonlinear system identification, Symmetry, Neural networks, Variational
integrators

††footnotetext: These authors contributed equally to this work.

It is a well-established theory from the field of variational calculus that
the behaviour of an unforced system can be derived from the
Lagrange-d’Alembert principle using only knowledge of the system’s Lagrangian.
The resulting Euler-Lagrange equations of motion act as a model describing how
the system’s configuration and velocity evolve as time progresses from an
initial displacement and thus serve as a state space system description.
Equivalently, one could describe the system in the phase space through the
canonical equations of motion based only on knowledge of the system’s
Hamiltonian. To learn a system representation, one could identify the
equations of motion directly, or alternatively, identify the Hamiltonian or
the Lagrangian which define them. Using the latter approach Greydanus et al.
(2019) proposed the Hamiltonian Neural Network (HNN), closely followed by the
Lagrangian Neural Networks (LNN) by Cranmer et al. (2020) leading to a rapid
growth of research interest in this topic.

Drawing from the field of mathematical physics we know that Hamiltonian and
Lagrangian models encode additional information about the system such as
symmetries. Through Noether’s theorem, these symmetries relate to conservation
laws which further restrict the motion of the system. The use of specialised
integrators, known as geometric integrators, allows the preservation of these
characteristics from the continuous into the discrete domain. This results in
a better qualitative representation of the system and a preservation of the
symmetries and energy accurately or up to small bounded oscillations for
exponentially long times (Hairer et al., 2006).

Variational integrators are geometric integrators that are based on the theory
of discrete mechanics (Marsden and West, 2001). With standard integrators, the
equations of motion are first derived and then discretised. With variational
integrators the discretisation is done one step earlier. In this case, a
discrete approximation of the Lagrangian is created and then used in a
discrete Lagrange-d’Alembert principle leading to equations of motion which
are structure-preserving. With this in mind Saemundsson et al. (2020) proposed
learning of the Lagrangian with forward simulations using the variational
integrator and demonstrated accurate long-term predictions even when learning
from noisy data. This was extended for systems with external forcing by Havens
and Chowdhary (2021). Both of these works rely on a phase space formulation of
the variational integrator, which requires data and prediction of both the
configuration and the momentum.

However, when using a variational integrator, we can use a second formulation
which removes the need to store observations or calculate predictions of the
canonical momentum or the velocity in both the learning and testing stages. In
this formulation an initial momentum value and the observations and
predictions of the configuration are sufficient to model the system behaviour
in time. This is an advantage compared to the use of standard integrators as
noted by Aoshima et al. (2021) for an alternative energy-preserving method in
the absence of external forcing. Another important competing approach are the
symplectic recurrent neural networks (SRNN) by Chen et al. (2019), which learn
the Hamiltonian and rely on a different symplectic integrator that is
performant even in noisy systems. However, it again relies on observations
from both the configuration and the momentum. The approaches with variational
integrators mentioned so far learn a continuous Lagrangian expression from
discrete data and use a discrete approximation of this expression in the
forward variational simulations. These approximations introduce discretisation
errors in the results which could amplify the modelling error of a learnt
Lagrangian (Ober-Blöbaum and Offen (2023)). To prevent these errors one could
instead directly learn the discrete Lagrangian and use it for variational
simulations in testing and training as demonstrated by Qin (2020) for
modelling nonlinear oscillations and the Kepler problem and by Santos et al.
(2022) for an assumed structure of a potential and a kinetic energy term
learned with two separate neural networks.

Despite these investigations, the question still remains: once you learn
either the continuous or the discrete Lagrangian how can one best obtain one
from the other? This question is non-trivial and goes back to the field of
backward error analysis (Hairer et al., 2006). When using a symplectic
integrator, one does not preserve the energy exactly, but rather preserves a
modified energy function. Thus, the motion of this predicted discrete
trajectory matches the continuous dynamics governed by a modified Lagrangian
instead of the true Lagrangian (see Figure 1). Offen and Ober-Blöbaum
investigated how to obtain such Hamiltonian and Lagrangian functions which
when discretised by the integrator would lead to the exact motion of the true
system model being predicted and the true conservation law being preserved.
They achieved this using the principle of variational backward error analysis
(VBEA) (Vermeeren, 2017) and named these functions inverse modified
Lagrangians (Ober-Blöbaum and Offen (2023); Offen and Ober-Blöbaum (2022)).
Their models were built using Gaussian Processes and allowed for Lagrangians
of any form, without assuming separable terms for potential and kinetic
energy. VBEA was also used to improve the accuracy in computing the initial
conjugate momentum and an additional normalisation was introduced into the
scheme which would prevent techniques such as the LNN from learning constant
Lagrangians; which satisfy the equations of motion but do not produce any
dynamics. Together these improvements allowed for the use of larger step sizes
in simulations while retaining high accuracy in the solution.

In this article we redefine their work and create a neural network model of
the discrete Lagrangian which consolidates and further improves upon the
advantages of previous approaches. Our approach aims to learn an inverse
modified discrete Lagrangian without prior restriction on its form. As
previously discussed the use of discrete inverse modified Lagrangians in
forward simulations results in a discrete prediction of the exact dynamics
rather than modified ones. Its discrete form avoids discretisation error and
requires observations and predictions of the configuration alone, without the
need to obtain or simulate the velocity as well. Furthermore, we propose an
additional degeneracy cost term which inductively biases the network to avoid
learning discrete Lagrangians whose equations of motion have degenerate roots
and pose problems for forward simulation incorporating numerical root-solving
techniques. The additional term also inductively biases the neural network to
avoid learning non-desirable constant Lagrangian solutions. With these
improvements, the first of our two proposed methods, the Discrete Lagrangian
Neural Network (DLNN) successfully incorporates the use of variational
integrators and variational backward error analysis in a novel neural network
based model, which uses the underlying structure of the problem to improve the
representation qualitatively and explicitly guard against numerical issues in
forward integration.

Inductive biases borne out of the incorporation of symmetries into learning
algorithms have previously shown desirable behaviours in the resulting models
such as reduced sample complexity and improved generalisation whilst
significantly reducing model complexity (Dehmamy et al., 2021). We introduce a
novel method to automatically learn variational symmetries of a dynamical
system along with the system’s discrete Lagrangian. More precisely, we
introduce a framework that for a given Lie group action on the configuration
space identifies a subgroup which acts by symmetries. Through construction of
a momentum map, conserved quantities are derived. The framework is explicitly
worked out for the group of affine linear transformations, which consists of
arbitrary compositions of translation, rotations, and scaling transformations.
We demonstrate in numerical examples that our proposed Symmetric Discrete
Lagrangian Neural Network (SymDLNN) successfully learns symmetries and
conservation laws and that incorporation of symmetries improves the
predictions compared to previous approaches.

To summarise, our contributions are:

  * •

A new regularisation term for learning discrete Lagrangians from snapshot data
of motions. The regularisation guarantees that the learnt discrete Lagrangian
can be efficiently employed in the numerical computation of trajectories.

  * •

A novel framework to automatically identify symmetries and conservation laws
while a discrete Lagrangian is learned.

  * •

The use of variational backward error analysis for learning discrete
Lagrangians and symmetries.

##  1 Background

###  1.1 Continuous Lagrangian dynamics

Consider a mechanical system, whose configuration $q(t)\in\mathbb{R}^{n_{q}}$
evolves on a configuration manifold $Q$, with associated tangent bundle $TQ$.
The behaviour of the system in time can be described in state space by the
evolution of its configuration vector $\underline{q}(t)\in Q$ and its
associated velocity vector $\underline{\dot{q}}(t)$ such that
$(\underline{q}(t),\underline{\dot{q}}(t))\in TQ$. The system is assumed to
possess a regular Lagrangian $\mathcal{L}:TQ\to\mathbb{R}$, which is not
explicitly time dependent. It is known that the motion of the system is
governed by the Lagrange-d’Alembert principle, which requires that

$\delta\int_{t_{0}}^{t_{f}}\mathcal{L}(\underline{q},\underline{\dot{q}})\,dt=0$

is satisfied for all variations $\delta q$ with $\delta q(t_{0})=\delta
q(t_{f})=0$ (Marsden and West, 2001). Using integration by parts this
principle results in the following equations:

$D_{1}\mathcal{L}(\underline{q},\underline{\dot{q}})-\frac{d}{dt}D_{2}\mathcal{L}(\underline{q},\underline{\dot{q}})=0$
(1)

known as the Euler-Lagrange equations of motion (Liberzon, 2011). Here $D_{i}$
is the partial derivative operator with respect to the $i$-th argument and the
dependencies on time have been omitted to simplify the notation. Using these
equations one can model how the system evolves in the state space over time.
The system’s evolution can equivalently be described in phase space using the
configuration vector $\underline{q}$ and conjugate momenta $\underline{p}$
defined as
$\underline{p}=D_{2}\mathcal{L}(\underline{q},\underline{\dot{q}})$. Notably,
that for a given system, the Lagrangian, which satisfies Equation (1) is not
unique and thus the conjugate momentum is not unique either. For example, any
constant function could satisfy the equation but would not lead to any
dynamics. This is important for schemes attempting to learn the Lagrangian
function as one would need to safeguard against constant solutions.

###  1.2 Discrete Lagrangian mechanics

In the discrete setting $TQ$ is replaced with $Q\times Q$ and a discrete
configuration path is defined as
$\underline{q}_{\,d}(t_{k})=\underline{q}_{\,k}$ with
$\underline{q}_{k}\approx\underline{q}(t_{k})$ for $t_{k}=t_{0}+k\Delta
T\textrm{ and }k=0,\ldots,N$ where $N=t_{f}/\Delta t$. Based on this
discretisation a discrete Lagrangian $\mathcal{L}_{d}$ can be defined as

$\mathcal{L}_{d}(\underline{q}_{k},\underline{q}_{k+1})\approx\int_{t_{k}}^{t_{k+1}}\mathcal{L}(\underline{q},\underline{\dot{q}})\,dt$

and plugged in a discrete version of the Lagrange-d’ Alembert principle:

$\delta\sum_{i=0}^{N-1}\mathcal{L}_{d}(\underline{q}_{k},\underline{q}_{k+1})=0$

leading to the discrete Euler-Lagrange equations:

$D_{2}\mathcal{L}_{d}(\underline{q}_{k-1},\underline{q}_{k})+D_{1}\mathcal{L}_{d}(\underline{q}_{k},\underline{q}_{k+1})=0$
(2)

To parallel the continuous theory one can define the discrete conjugate
momentum $\underline{p}_{k}$ as
$\underline{p}_{k}=-D_{1}\mathcal{L}_{d}(\underline{q}_{k},\underline{q}_{k+1})$
(see Marsden and West (2001)). Given $\underline{q}_{0}$ and
$\underline{q}_{1}$ one can use equations (2) to time-step in time to compute
$\underline{q}_{k}$ for $k>1$. To compute $\underline{q}_{1}$ and initiate the
integration one can use the expression for the discrete momentum for the
initial position $\underline{q}_{0}$ and conjugate momentum
$\underline{p}_{0}$. Similarly to the continuous section techniques learning
the discrete Lagrangian would need to safeguard against solutions which would
satisfy Equation (2) but would not produce meaningful dynamics.

###  1.3 Standard and variational backward error analysis

Assume we have a system with continuous dynamics $\dot{x}(t)=f(x(t))$ and a
numerical integrator $\dot{x}_{k+1}=z(x_{k})$. This discrete method is only an
approximation of the actual dynamics $x(k\Delta t)$ at time $t=k\Delta t$ or
in other words $x_{k}\approx x(k\Delta t)$. Backward error analysis focuses on
discovering dynamics $\dot{x}^{m}=f^{m}(x^{m}(t))$ for which equation
$\dot{x}_{k+1}=z(x_{k})$ is in fact the exact solution or in other words
$x_{k}=x^{m}(k\Delta t)$. The dynamics $\dot{x}^{m}=f^{m}(x^{m}(t))$ are often
termed as the modified equations. By studying the difference between $x_{k}$
and $x^{m}(k\Delta t)$ one can determine important properties of the
integrator scheme (Hairer et al., 2006).

Following from this idea, Vermeeren (2017) has developed a similar analysis
for variational integrators, known as variational backward error analysis. As
depicted in Figure 1, the true dynamics $\underline{q}(t)$ of the system are
governed by the continuous equations (1) when the true system Lagrangian is
used $\mathcal{L}^{true}$. The solutions of Equation (2)
$\\{\underline{q}^{m}_{k}\\}_{k=0}^{N}$ are only approximations of these
dynamics. On the other hand Vermeeren (2017) derived a modified Lagrangian
expression $\mathcal{L}^{mod}$ which when plugged in (1) would obtain a
solution $\underline{q}^{m}(t)$ such that $\underline{q}^{m}(k\Delta
t)=\underline{q}^{m}_{k}$. Ober-Blöbaum and Offen (2023); Offen and Ober-
Blöbaum (2022) on the other hand searched for an inverse modified Lagrangian
$\mathcal{L}^{invmod}$ expression which, when plugged into the variational
scheme would obtain discrete path $q_{k}$ which accurately represents the true
dynamics in the discrete domain i.e. $\underline{q}_{k}=\underline{q}(k\Delta
t)$. He further describes how $\mathcal{L}^{invmod}$ can be used to obtain
$\mathcal{L}^{true}$.

Figure 1: Diagramatic explanation of the modified and inverse modified
Lagrangian functions (contin.-continuous, EOM- equations of motion).

In this article, we propose to learn discrete inverse modified Lagrangians
directly from data. Once it is learned, motions can be computed using Equation
(2). To initialise the computation, initial data consisting of the first two
positions $(q_{0},q_{1})$ of the trajectory are required. For the purpose of
system identification, the inverse modified Lagrangian $\mathcal{L}^{invmod}$
can be computed from its discrete learnt counterpart by a change of variables
depending on the variational integrator and subsequently using variational
backward error analysis formulas used to identify a function for the true
Lagrangian, which will be denoted $\mathcal{L}^{VBEA}$.

For example, for a simulation using a step of size $\Delta t$ and the
variational midpoint rule in the one-dimensional case
$\mathcal{L}^{invmod}(\underline{q},\underline{\dot{q}})=\mathcal{L}_{d}^{invmod}(\underline{q}-\Delta
t/2\dot{\underline{q}},\underline{q}+\Delta t/2\dot{\underline{q}})$, where
$\mathcal{L}_{d}^{invmod}$ is the discrete Lagrangian that we learn, and

$\displaystyle\mathcal{L}^{VBEA}$
$\displaystyle=\mathcal{L}^{invmod}+\tfrac{1}{24}\Delta
t^{2}\frac{(\partial_{q}{\mathcal{L}^{invmod}}-\partial^{2}_{q\dot{q}}{\mathcal{L}^{invmod}}\dot{q})^{2}}{\partial^{2}_{\dot{q}\dot{q}}{\mathcal{L}^{invmod}}}$
$\displaystyle-\tfrac{1}{24}\Delta
t^{2}\partial_{qq}{\mathcal{L}^{invmod}}\dot{\underline{q}}^{2}+O(\Delta
t^{4}).$ (3)

Formulas for higher dimensions or truncation orders can be derived as
explained in Vermeeren (2017).

###  1.4 Variational symmetries, conserved quantities, and a framework for our
symmetry learning method

Dynamical systems governed by Euler-Lagrange equations can exhibit symmetries.
These are of great significance since they encode important qualitative
features of motions: for instance, the idealised motion of a pendulum on a
cart is described by Equation (1) with the Lagrangian

$\mathcal{L}_{\mathrm{cp}}(s,\phi,\dot{s},\dot{\phi})=\frac{1}{2}(\alpha\dot{\phi}^{2}+2\beta\cos(\phi)\dot{s}\dot{\phi}+\gamma\dot{s}^{2})+D\cos(\phi),$
(4)

with position variables $q=(s,\phi)$ and velocity variables
$\dot{q}=(\dot{s},\dot{\phi})$, where $\alpha=m_{1}l^{2}$, $\beta=m_{1}l$,
$\gamma=m_{2}+m_{1}$, and $D=-m_{1}gl$. Here $m_{2}$ is the mass of the cart,
$m_{1}$ the mass of the pendulum, $l$ the length of the pendulum, $g$ a
gravity acceleration constant. The Lagrangian $\mathcal{L}_{\mathrm{cp}}$ does
not depend explicitly on the variable $s$. In other words, it is invariant
under translations of $s$, that is, the action $s\mapsto s+\eta$ for
$\eta\in\mathbb{R}$ and $\frac{\partial\mathcal{L}_{\mathrm{cp}}}{\partial
s}(\underline{q},\underline{\dot{q}})=0$ for all
$(\underline{q},\underline{\dot{q}})$. By Noether’s theorem, this symmetry
corresponds to the conservation of the conjugate momentum

$I_{\mathrm{cp}}=\frac{\partial{\mathcal{L}}_{\mathrm{cp}}}{\partial\dot{s}}=\beta\cos(\phi)\dot{\phi}+\gamma\dot{s}$
(5)

along motions $q(t)=(s(t),\phi(t))$. More generally, if for a Lagrangian
$\mathcal{L}(\underline{q},\underline{\dot{q}})$ with $\underline{q}\in
Q=\mathbb{R}^{n_{d}}$ there exists a direction $w\in\mathbb{R}^{n_{d}}$ such
that $\mathcal{L}$ is invariant under an action
$(\underline{q},\underline{\dot{q}})\mapsto(\underline{q}+sw,\dot{\underline{q}})$
for all $s\in\mathbb{R}$ then the directional derivative

$w^{\top}\nabla_{\underline{q}}\mathcal{L}(\underline{q},\underline{\dot{q}})=\sum_{j=1}^{n_{d}}w_{j}\frac{\partial\mathcal{L}}{\partial
q_{j}}(\underline{q},\underline{\dot{q}})=0$ (6)

and the quantity

$I(\underline{q},\underline{\dot{q}})=w^{\top}\nabla_{\dot{q}}\mathcal{L}(\underline{q},\underline{\dot{q}})=\sum_{j=1}^{n_{d}}w_{j}\frac{\partial{\mathcal{L}}}{\partial\dot{q}_{j}}(\underline{q},\underline{\dot{q}})$
(7)

are preserved along motions $\underline{q}(t)$ of Equation (1). Even more
generally, an affine linear transformation is described by an invertible
matrix $W\in\mathrm{Gl}(\mathbb{R},n_{q})$ and a vector
$w\in\mathbb{R}^{n_{q}}$ and acts on $(\underline{q},\underline{\dot{q}})\in
TQ$ by
$(\underline{q},\underline{\dot{q}})\mapsto(W\underline{q}+w,W\underline{\dot{q}})$.
Here $\mathrm{Gl}(\mathbb{R},n_{q})$ denotes the group of invertible matrices.
The matrix $W$ encodes a rotation and scaling, and the vector $w$ encodes a
translation. We can write $W$ and $w$ into a matrix $\begin{pmatrix}W&w\\\
0&1\end{pmatrix}\in\mathrm{Gl}(\mathbb{R},n_{q}+1)$. Affine transformations
form a group $G$ which can be represented by the following subgroup of
$(n_{q}+1)\times(n_{q}+1)$-dimensional invertible matrices

$G_{\mathrm{aff}}=\left\\{\left.\begin{pmatrix}W&w\\\
0&1\end{pmatrix}\,\right|\,A\in\mathrm{Gl}(\mathbb{R},n_{q}),\,w\in\mathbb{R}^{n_{q}}\right\\}.$

Any 1-dimensional subgroup of $G_{\mathrm{aff}}$ can be defined by a matrix
$M\in\mathbb{R}^{n_{q}\times n_{q}}$ and vector $w\in\mathbb{R}^{n_{q}}$ and
is of the form

$G_{(M,w)}=\left\\{\left.\exp\left(\begin{pmatrix}\eta M&\eta w\\\
0&1\end{pmatrix}\right)\right|\eta\in\mathbb{R}\right\\}\subset G,$

where $\exp$ is the matrix exponential. If $G_{(M,w)}$ acts by symmetries,
that is, if the Lagrangian $\mathcal{L}$ is invariant under actions by the
elements in $G_{(M,w)}$, then

$(M\underline{q}+w)^{\top}\nabla_{\underline{q}}\mathcal{L}(\underline{q},\underline{\dot{q}})+M^{\top}\nabla_{\underline{\dot{q}}}\mathcal{L}(\underline{q},\underline{\dot{q}})=0$
(8)

for all $(\underline{q},\underline{\dot{q}})\in TQ$ and the quantity

$I(\underline{q},\underline{\dot{q}})=(M\underline{q}+w)^{\top}\nabla_{\underline{\dot{q}}}\mathcal{L}(\underline{q},\underline{\dot{q}})$
(9)

is conserved along motions defined by the Euler-Lagrange Equation (1). When
working with discrete Lagrangians $\mathcal{L}_{d}\colon Q\times
Q\to\mathbb{R}$, the symmetry condition (8) is replaced by

$(M{\underline{q}}_{0}+w)^{\top}\nabla_{\underline{q}_{0}}\mathcal{L}({\underline{q}}_{0},{\underline{q}}_{1})+(M{\underline{q}}_{1}+w)^{\top}\nabla_{\underline{q}_{1}}\mathcal{L}({\underline{q}}_{0},{\underline{q}}_{1})=0$
(10)

for all $(\underline{q}_{0},\underline{q}_{1})\in Q\times Q$ and the conserved
quantity for discrete motions becomes

$I({\underline{q}}_{k},{\underline{q}}_{k+1})=-(M{\underline{q}}_{k}+w)^{\top}\nabla_{q_{k}}\mathcal{L}({\underline{q}}_{k},{\underline{q}}_{k+1}).$
(11)

For details, we refer to the book by Marsden and Ratiu (1999). Our method
SymDLNN automatically identifies subgroups $G_{(M,w)}$ of the group of affine
linear transformations under which a (discrete) Lagrangian is invariant while
the (discrete) Lagrangian of a system is learned. As a consequence, we
identify the corresponding conserved quantity as well.

Lagrangians are not uniquely determined by the system’s motion and non-
symmetric Lagrangians can govern highly symmetric dynamical systems. Our
method will guide the learning process to a symmetric Lagrangian. This
regularises the learning process and improves the reliability of predictions
as we will demonstrate in numerical examples. Our approach is not restricted
to affine linear symmetries or 1-dimensional symmetry groups. To present the
most general framework, let us briefly introduce Lie group actions and
invariant vector fields. See Marsden and Ratiu (1999) for details. For a Lie
group $G$, let $\mathfrak{g}$ denote its Lie algebra and
$\exp\colon\mathfrak{g}\to G$ the exponential map. Consider a group action
$a\colon G\to\mathrm{Diff}(Q),\,g\mapsto a_{g}$. Here $\mathrm{Diff}(Q)$
denotes the group of diffeomorphisms on $Q=\mathbb{R}^{n_{q}}$. In the setting
of continuous Lagrangians, the group action can be prolonged to an action
$A\colon G\to\mathrm{Diff}(TQ)$, $g\mapsto A_{g}$ by defining
$A_{g}(\underline{q},\underline{\dot{q}})=(a_{g}({\underline{q}}),Da_{g}({\underline{q}})\dot{{\underline{q}}})$,
where $Da_{g}({\underline{q}})$ is the Jacobi matrix of the diffeomorphism
$a_{g}$ at $q$. Here we have identified
$\mathcal{M}=TQ\cong\mathbb{R}^{n_{q}}\times\mathbb{R}^{n_{q}}$. For the
setting of discrete Lagrangians, the diagonal group action $A\colon
G\to\mathrm{Diff}(Q\times Q)$, $g\mapsto A_{g}$ with
$A_{g}(q_{0},q_{1})=(a_{g}(q_{0}),a_{g}(q_{1}))$ is considered instead of the
prolonged group action and $\mathcal{M}=Q\times Q$. In both cases, for
$v\in\mathfrak{g}$ the left invariant vector field
$\hat{v}\in\mathfrak{X}(\mathcal{M})$ is defined by

$\hat{v}_{z}=\left.\frac{\d{missing}}{\d{t}}\right|_{t=0}A_{\exp{(tv)}}(z)\in
T_{z}\mathcal{M},\quad z\in\mathcal{M}.$

These vector fields can be thought of as infinitesimal actions of the Lie
group $G$ on $\mathcal{M}$.

The idea of SymDLNN is to identify a basis $v^{1},\ldots,v^{K}$ of a subspace
of $V^{(K)}\subset\mathfrak{g}$ such that $\hat{v}^{j}_{z}(\mathcal{L})=0$ for
all $z\in\mathcal{M}$ and $j=1,\ldots,K$. Under mild assumptions on the group
action, a momentum map $J\colon M\to V^{(K)\ast}$ can be constructed from
which $K$ functionally independent conserved quantities can be computed as
$I^{j}(\underline{q},\underline{\dot{q}})=\langle
J(\underline{q},\underline{\dot{q}}),v^{j}\rangle$, where
$\langle\cdot,\cdot\rangle$ denotes the dual pairing. Relating this theory
back to the example of affine linear symmetries, the Lie algebra for the
affine linear group $G_{\mathrm{aff}}$ can be written as the following
subspace of $\mathbb{R}^{(n_{q}+1)\times(n_{q}+1)}$:

$\mathfrak{g}=\left\\{\left.\begin{pmatrix}M&w\\\
0&0\end{pmatrix}\,\right|\,M\in\mathbb{R}^{n_{q}\times
n_{q}},\,w\in\mathbb{R}^{n_{q}}\right\\}$ (12)

SymDLNN seeks as many elements $v^{j}=\begin{pmatrix}M^{j}&w^{j}\\\
0&0\end{pmatrix}\in\mathfrak{g}$ as possible for which the symmetry condition
(8) (or (10) in case of discrete Lagrangians) holds and such that all
$v^{1},\ldots,v^{(K)}$ are linearly independent. Then (under non-triviality
conditions on $\mathcal{L}$) the $K$ quantities provided by Equation (9) or
Equation (11) are functionally independent and are conserved under motions.

##  2 Proposed approach

Our scheme aims to obtain a model for the true Lagrangian of the system with
no restriction on the structure of the Lagrangian function apart from its
regularity. This is done by learning a discrete Lagrangian function
$\mathcal{\mathcal{L}}^{NN}_{d}$ which corresponds to the discrete inverse
modified Lagrangian $\mathcal{L}^{invmod}_{d}$. After the learning stage is
completed, these functions are used for forward simulation with a variational
integrator to obtain an accurate discretisation $\underline{q}_{k}$ of the
true dynamics. The continuous inverse modified Lagrangian is obtained using
$\mathcal{L}^{invmod}(\underline{q},\underline{\dot{q}})=\mathcal{\mathcal{L}}^{NN}_{d}(\underline{q}-\Delta
t\underline{\dot{q}}/2,\underline{q}+\Delta t\underline{\dot{q}}/2)$ based on
the variational midpoint rule and employing the theory of VBEA from Section
1.3, an expression of the true Lagrangian $\mathcal{L}^{VBEA}$ is obtained and
used to identify the corresponding Hamiltonian $\mathcal{H}^{VBEA}$.

To learn $\mathcal{L}^{NN}_{d}$ we rely on observations and simulations only
of the configuration in time and do not pose any restrictions on its form. If
required, velocity data at step $k$ can be computed using central difference
or by solving

$-D_{1}\mathcal{L}^{NN}_{d}(\underline{q}_{k},\underline{q}_{k+1})=\nabla_{\dot{q}}\mathcal{L}^{VBEA}(\underline{q},\underline{\dot{q}})|_{(}\underline{q}_{k},\underline{\dot{q}}_{k})$

for $\underline{\dot{q}}_{k}$. Observations and predictions are created from
configuration data using three consecutive values in discrete time
$(\underline{q}_{k-1},\underline{q}_{k},\underline{q}_{k+1})$ from the same
test trajectory. Each triple is used to compute
$D_{2}\mathcal{L}_{d}(\underline{q}_{k-1},\underline{q}_{k})+D_{1}\mathcal{L}_{d}(\underline{q}_{k},\underline{q}_{k+1})$
As known from Equation (1) these values should be zero for all discrete points
of the trajectory. Thus the neural network can be trained by minimising

$\min\frac{1}{N}\sum_{k=0}^{N-1}\sum_{i=1}^{n_{q}}(D_{2}\mathcal{L}_{d}(\underline{q}_{k-1},\underline{q}_{k})+D_{1}\mathcal{L}_{d}(\underline{q}_{k},\underline{q}_{k+1}))[i])^{2}$
(13)

with a gradient descent algorithm (here $[\cdot]_{i}$ denotes the $i$-th
element of the vector).

However, using only Equation (13) without properly safeguarding the learning
process can learn a model in which the roots of the equations of motion (the
trajectory values) are degenerate, only touching the zero axis and not
necessarily crossing it. Such roots present difficulties for forward
simulation as they are difficult to find using common numerical methods such
as Newton-Rhapson. To prevent such behaviour we propose adding the degeneracy
term

$\displaystyle\sum_{k=0}^{N-1}\Big{(}1-\frac{1}{1-e^{-0.01d_{k}}}\Big{)},\;\;d_{k}^{m}=\textrm{det}\Big{(}D_{2}D_{1}\mathcal{L}_{d}(\underline{q}_{k},\underline{q}_{k+1})\Big{)}$

to the cost function (13), where $m=1\textrm{ or }2$ and $d$ corresponds to
the Jacobian of the discrete equations of motion (2). This term aims to
increase the slope of the equations of motion at the root locations,
safeguarding against degenerate roots (which have vanishing gradients at the
root location) and rewarding steeper crossings. This term also prevents the
learning of constant discrete Lagrangian functions, which satisfy the
equations of motion but do not result in meaningful dynamics.

###  2.1 Incorporation of symmetry into the loss function

Using the proposed loss function we learn a non-degenerate discrete inverse
modified Lagrangian, which can then be used to predict motions. However, since
discrete Lagrangians are not uniquely determined by the system’s motion, the
learned discrete Lagrangian for a dynamical system which exhibits symmetries
can be arbitrarily unsymmetric even if the learnt discrete Lagrangian
minimises the loss function. In the following, we introduce a method which
automatically detects symmetries of the dynamical system and drives the learnt
discrete Lagrangian towards a symmetric representation. This is useful for
system identification as symmetries inform us about conserved quantities which
constitute important qualitative features of the dynamical system. Moreover,
driving the learnt Lagrangian to a symmetric representation acts as an
additional regulariser in the learning process and can be beneficial for
numerical simulations of the learned system.

Based on the symmetry condition (10) of Section 1.4, to detect an affine
linear symmetry of the dynamical system we propose adding the following term
to the cost function:

$\displaystyle\ell_{\mathrm{sym}}=\frac{1}{N}\sum_{k=0}^{N-1}|(M\underline{q}_{k}w)^{\top}\nabla_{\underline{q}_{k}}\mathcal{L}_{d}^{NN}(\underline{q}_{k},\underline{q}_{k+1})$
$\displaystyle+(M\underline{q}_{k+1}+w)^{\top}\nabla_{\underline{q}_{k+1}}\mathcal{L}_{d}^{NN}(\underline{q}_{k},\underline{q}_{k+1})|^{2},$

where $\underline{q}_{0},\ldots,\underline{q}_{N}$ is a training trajectory.
Additionally, the above expression is summed over all training trajectories.
The $n_{q}^{2}$ elements in the matrix $M$ and the $n_{q}$ elements in $w$ are
now trainable parameters in addition to the parameters of
$\mathcal{L}_{d}^{NN}$ and are included in the minimisation process.
Furthermore, we add the non-triviality condition

$|\|M\|^{2}+\|w\|^{2}-1|^{2}$ (14)

to the loss function to encourage the process to learn a non-trivial symmetry.
After learning, Equation (11) is a candidate for a conserved quantity of the
discrete system and Equation (9) of the underlying continuous system.

###### Remark 1

The choice to test the symmetry condition (10) on points
$({\underline{q}}_{k},{\underline{q}}_{k+1})$ in the training data to
construct $\ell_{\mathrm{sym}}$ is arbitrary. In principle,
$\ell_{\mathrm{sym}}$ could be any approximation (e.g. a Monte-Carlo
integration) of

$\displaystyle\frac{1}{\mathrm{vol}(\mathcal{M}^{\mathrm{o}})}\int_{\mathcal{M}^{\mathrm{o}}}|(M{\underline{q}}_{0}+w)^{\top}\nabla_{\underline{q}_{0}}\mathcal{L}_{d}^{NN}({\underline{q}}_{0},{\underline{q}}_{1})$
$\displaystyle+(M{\underline{q}}_{1}+w)^{\top}\nabla_{q_{1}}\mathcal{L}_{d}^{NN}({\underline{q}}_{0},{\underline{q}}_{1})|^{2}\d{\text@underline{q}}_{0}\d{\text@underline{q}}_{1}$

where $\mathcal{M}^{\mathrm{o}}\subset Q\times Q$ is a (topologically open,
pre-compact) subset of the discrete phase space covering all parts of
interest.

###### Remark 2

To learn $K$ functionally independent integrals of motions, we simply add $K$
instances $\ell_{\mathrm{sym}}^{(j)}$ of $\ell_{\mathrm{sym}}$ from Equation
(15) to the loss function. Each $\ell_{\mathrm{sym}}^{(j)}$ has trainable
parameters $(M^{(j)},w^{(j)})$. Further, the non-triviality condition (14) is
added for each instance of $\ell_{\mathrm{sym}}^{(j)}$. To make sure that all
learned $(M^{(k)},w^{(k)})$ yield functionally independent integrals of
motions, we need to make sure that they span a basis of the Lie algebra (12).
This is achieved by adding the orthogonality condition

$\sum_{k=1}^{K}\sum_{s=1}^{k-1}(\mathrm{vec}(M^{(s)})^{\top}\mathrm{vec}(M^{(k)})+(w^{(s)})^{\top}w^{(k)})$

to the loss function. Here $\mathrm{vec}$ writes the columns of a matrix into
a single vector for the computation of the Frobenius inner product of two
matrices. This corresponds to learning a $K$-dimensional symmetry group.

###### Remark 3

The approach can be generalised to arbitrary Lie group actions as considered
at the end of Section 1.4 as follows. For $k=1,\ldots,K$ define
$\ell_{\mathrm{sym}}^{(k)}$ as a numerical approximation of

$\ell_{\mathrm{sym}}^{(k)}\approx\frac{1}{\d{\mathrm{missing}}{vol}(\mathcal{M}^{\mathrm{o}})}\int_{\mathcal{M}^{\mathrm{o}}}|\hat{v}^{(k)}(\mathcal{L}_{d})|^{2}\d{\text@underline{q}}_{0}\d{\text@underline{q}}_{1}$
(15)

Here $\hat{v}^{(k)}$ denotes the invariant vector field to $v$. The term
$\ell_{\mathrm{sym}}^{(k)}$ measures how invariant $\mathcal{L}$ is under
actions with group elements of $\exp({tv^{(k)}}|t\in\mathbb{R})$. Equip
$\mathfrak{g}$ with an inner product $\langle\cdot,\cdot\rangle$ and norm
$\|\cdot\|$. Given weights $\alpha^{(k)},\beta^{(k)}>0$ define

$\displaystyle\ell^{\mathrm{total}}_{\mathrm{sym}}=\sum_{k=1}^{K}\left(\ell_{\mathrm{sym}}^{(k)}+\alpha^{(k)}|\|v^{(k)}\|-1|^{2}+\beta^{(k)}\sum_{s=1}^{k-1}\langle
v^{(k)},v^{(s)}\rangle\right)$

which is added to the loss function. Here $\alpha^{(k)}$ and $\beta^{(k)}$ are
fixed, non-negative weights. The last two terms of $\ell_{\mathrm{sym}}$
measure the orthonormality of the spanning set $v^{(1)},\ldots,v^{(K)}$ while
the first term measures how well infinitesimal actions by elements of
$\exp(V)$ preserve $\mathcal{L}$.

##  3 Results

To evaluate the proposed DLNN and SymDLNN methods alongside LNN, each of the
methods was applied to two systems: a pendulum on a cart defined by the
Lagrangian in Equation (4) with $m_{1}=m_{2}=l=1$, and the translational
symmetry with $M=[0,0;0,0]$, $w=[1,0]$ and the Kepler problem in two
dimensions with Lagrangian

$\mathcal{L}_{\mathrm{Kp}}(x,y,\dot{x},\dot{y})=\frac{1}{2}(\dot{x}^{2}+\dot{y}^{2})+\frac{Gm_{1}m_{2}}{\sqrt{x^{2}+y^{2}}}$
(16)

for $\underline{q}=[x,y]^{T}$, $G=6.673\times 10^{-26}$, $m_{1}=6\times
10^{24}$, $m_{2}=100$ and rotational symmetry of the form $w=[0,0]$,
$M=[0,\sqrt{2}/2;-\sqrt{2}/2,0]$. For each system, the methods are tasked with
learning the underlying dynamics defined by the Lagrangian in three evaluation
scenarios: a single trajectory case, a multi-trajectory case, and learning a
single trajectory with noisy measurements. After training a variational
integrator is used for simulations with each of the three methods in order to
assess how well the resulting system model is capable of: recreating the
trajectory it was trained on, its ability to predict and extension of this
trajectory and ability to preserve the energy and symmetry. Code for all
implementations and experiments can be found at
https://github.com/yanalish/SymDLNN.

Figure 2: Pendulum on the cart example: Recreation and prediction of the
trajectory based on single trajectory observation. Figure 3:  Kepler example:
Recreation and prediction of the trajectory based on single trajectory
observation.

In the first scenario a set of experiments are conducted wherein the
observations used to train the neural network are obtained from a single
trajectory of $N$ consecutive configuration points spaced at a stepsize of
$\Delta t$. For the pendulum on a cart $N=200$ with $\Delta t=0.01$ was
chosen, whereas for the Kepler problem these were set to $N=50$ and $\Delta
t=0.1$ Trajectory observations were obtained using a variational integrator
with a fine grained step size $\Delta t=10^{-4}$ for the pendulum on a cart
and $\Delta t=10^{-3}$ for the Kepler problem. Using these observations each
method attempts to recreate the trajectory with the learnt Lagrangian and
predict its extension for an additional $N_{extra}=100$ steps for the pendulum
and $N_{extra}=150$ for the Kepler problem, the resulting trajectory will be
denoted as $\underline{q}^{NN}$. For a fair comparison each of LNN, DLNN and
SymDLNN methods approximate the continuous and discrete Lagrangian
respectively using a 3-layer multi-layer perceptron (MLP) with 128 dimensional
hidden layers. The minimisation of the objective functions is performed
through gradient descent. Specifically, for all methods the optimisation is
performed using an Adam optimiser with an initial learning rate of $3\times
10^{-}3$, other Adam hyperparameters were set to $\beta_{1}=0.9$,
$\beta_{2}=0.999$, and $\epsilon=1\times 10^{-8}$. Each network was trained
for 100,000 epochs. The only notable difference in the training between the
methods is in the SymDLNN method which does not employ the symmetry term for
the first 5,000 epochs and then utilises them in the remaining 95,000 epochs.
The rationale for this is to initially optimise the parameters of the network
to find a suitable initial Lagrangian as the DLNN does and then discover the
symmetries.

In Figures 2 and 3 we can see for both example systems that the DLNN and
SymDLNN methods better recreate the original train trajectory than the LNN
approach. For the pendulum on a cart example DLNN provides the best prediction
whereas for Kepler both DLNN and SymDLNN outperform LNN for the unseen
extension. For the cart pendulum example the SymDLNN model identified symmetry
parameters $M^{NN}_{cp}=[1.51\times 10^{-3},-1.26\times 10^{-3};-7.83\times
10^{-6},-3.62\times 10^{-5}]$, $w^{NN}_{cp}=[9.97\times 10^{-1},8.04\times
10^{-2}]^{T}$ and for the Kepler $M^{NN}_{Kp}=[0.074,0.813;-0.568,-0.018]$,
$w^{NN}_{Kp}=[0.057,0.076]^{T}$ from initialization guess
$M^{guess}_{cp}=[0.1,0.1;0.1,0.1]$, $w^{guess}_{cp}=[1.5,0.5]^{T}$ and
$M^{guess}_{Kp}=[0.1,0.807;-0.607,0.1]$, $w^{guess}_{Kp}=[0.1,0.1]^{T}$
respectively. We believe this can be further improved with a more rigorous
hyperparameter search and more generous computing resources. Figures 4 and 5
demonstrate how accurately the symmetry and energy are preserved with each
model. In these plots, we compare how well LNN and DLNN preserve the symmetry
with

$\displaystyle
I^{true}_{k}=(\underline{p}_{k}^{NN})^{T}(M\underline{q}_{k}^{NN}+w)$ for
$\displaystyle\underline{p}_{k}^{NN}=\nabla_{\underline{\dot{q}}}\mathcal{L}^{true}\Big{(}\underline{q}_{k},\tfrac{\underline{q}_{k+1}^{NN}-\underline{q}_{k-1}^{NN}}{2\Delta
t}\Big{)}$

and how well SymDLLN learns it by computing

$\displaystyle
I^{NN}_{k}=(\underline{p}_{k}^{NN})^{T}(M^{NN}\underline{q}_{k}^{NN}+w^{NN})$
for
$\displaystyle\underline{p}_{k}^{NN}=-D_{1}\mathcal{L}^{NN}_{d}(\underline{q}^{NN}_{k},\underline{q}^{NN}_{k+1})$

For the energy plots we compare how well LNN has learnt the original
Hamiltonian and how accurately a Hamiltonian has been approximated based on
VBEA using the DLNN and SymDLNN model calculating respectively
$H^{true}_{k}(\underline{q}_{k},\tfrac{\underline{q}_{k+1}^{NN}-\underline{q}_{k-1}^{NN}}{2\Delta
t})$ and
$H^{VBEA}_{k}(\underline{q}_{k},\tfrac{\underline{q}_{k+1}^{NN}-\underline{q}_{k-1}^{NN}}{2\Delta
t})$. We can clearly see the improvement of using DLNN compared to the LNN
approach. The SymDLNN approach performs slightly worse, however, it is
important to note that for LNN and DLNN we plot based on our knowledge for the
symmetry expression whereas SymDLNN has learnt the symmetry itself, providing
a qualitatively better model and additional knowledge of the system’s
dynamics.

Figure 4: Pendulum on a cart example: Symmetry and energy error test based on
single trajectory observation (RP- recreation phase, PP- prediction phase).

Both examples represent non-chaotic systems, thus learning from a single
trajectory learns the behaviour of the system for a restricted portion of the
phase space. For this purpose, the same experiments were conducted for all
three methods when learning from data from 100 different trajectories of the
same lengths and steps as in the single trajectory examples. Recreation and
prediction of the trajectory as well as symmetry and energy preservation tests
closely resembled the results demonstrated above for learning from a single
trajectory.

Figure 5: Kepler example: Symmetry and energy error test based on single
trajectory observation (RP- recreation phase, PP- prediction phase). Figure 6:
Kepler example with noise: Recreation and prediction of the trajectory based
on observations gathered from a single trajectory. Figure 7: Kepler example
with noise: Symmetry and energy error test based on single trajectory
observation (RP- recregation phase, PP- prediction phase).

In our third and last evaluation scenario we examined the performance of our
approaches DLNN and SymDLNN in the presence of noise. As we can see below for
the Kepler problem the proposed model was capable of accurately recreating the
observed trajectory and extending it beyond for $N=50$ and $N_{extra}=150$
learning from a single trajectory observation for which measurement noise was
added with a standard Gaussian distribution with $\sigma^{2}=0.001$. For DLNN,
the $H^{VBEA}$ energy and the true symmetry are quite well preserved both
during the recreation and the prediction of the trajectory. Despite the noise,
SymDLNN is capable of preserving energy and although the symmetry oscillation
is slightly larger, it is based on a symmetry which the SymDLNN approach
itself identified with values
$M^{NN}=[0.047,0.814;-0.569,-0.005],w^{NN}=[0.064,0.079]^{T}$ with the same
initialization as before $M_{Kp}^{guess}$, $w_{Kp}^{guess}$. Moreover, we can
see that the additional symmetry term helps SymDLNN better extend the
trajectory than DLNN.

##  4 Conclusion

We have presented two novel methods, DLNN and SymDLNN, which successfully
incorporate the use of variational integrators and variational backward error
analysis for learning discrete Lagrangians, utilising the underlying structure
of the problem to improve the learned system representation qualitatively and
guard against numerical issues in forward simulations. This was achieved
through the proposal of a novel regularisation degeneracy term that
inductively biases the network to avoid learning Lagrangians whose equations
of motion have degenerate roots as well as constant Lagrangians. SymDLNN
further extends the DLNN in its ability to automatically identify symmetries
during the learning process through a framework that identifies a subgroup
which acts by symmetries upon a given Lie group on the configuration space.
Numerical examples on two example systems were used to demonstrate the
qualitative improvements of our proposed methods to previous approaches.

Y.L. acknowledges funding from the EPSRC Doctoral Training Partnership
EP/R513295/1, project reference 2280382. P.S. acknowledges funding by the W.D.
Armstrong Fund from the School of Technology at the University of Cambridge.
S.R. acknowledges funding by the EPSRC Centre for Doctoral Training in
Autonomous Intelligent Machines & Systems EP/L015897/1. C.O. acknowledges
funding by the Ministry of Culture and Science of the State of North Rhine-
Westphalia.

## References

  * Aoshima et al. (2021) Aoshima, T., Matsubara, T., and Yaguchi, T. (2021).  Deep discrete-time lagrangian mechanics.  In _ICLR2021 Workshop on Deep Learning for Simulation (SimDL)_ , volume 5.
  * Chen et al. (2019) Chen, Z., Zhang, J., Arjovsky, M., and Bottou, L. (2019).  Symplectic recurrent neural networks.  _arXiv preprint arXiv:1909.13334_.
  * Cranmer et al. (2020) Cranmer, M., Greydanus, S., Hoyer, S., Battaglia, P., Spergel, D., and Ho, S. (2020).  Lagrangian neural networks.  _arXiv preprint arXiv:2003.04630_.
  * Dehmamy et al. (2021) Dehmamy, N., Walters, R., Liu, Y., Wang, D., and Yu, R. (2021).  Automatic symmetry discovery with lie algebra convolutional network.  _Advances in Neural Information Processing Systems_ , 34, 2503–2515.
  * Greydanus et al. (2019) Greydanus, S., Dzamba, M., and Yosinski, J. (2019).  Hamiltonian neural networks.  _Advances in neural information processing systems_ , 32.
  * Hairer et al. (2006) Hairer, E., Lubich, C., and Wanner, G. (2006).  _Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations_ , volume 31.  vol. 31. Springer Science & Business Media, 2 edition.
  * Havens and Chowdhary (2021) Havens, A. and Chowdhary, G. (2021).  Forced variational integrator networks for prediction and control of mechanical systems.  In _Learning for Dynamics and Control_ , 1142–1153. PMLR.
  * Liberzon (2011) Liberzon, D. (2011).  _Calculus of variations and optimal control theory: a concise introduction_.  Princeton university press.
  * Marsden and Ratiu (1999) Marsden, J.E. and Ratiu, T.S. (1999).  _Introduction to mechanics and symmetry_ , volume 17 of _Texts in Applied Mathematics_.  Springer, 2nd edition.
  * Marsden and West (2001) Marsden, J.E. and West, M. (2001).  Discrete mechanics and variational integrators.  _Acta Numerica_ , 10(1), 357–514.
  * Ober-Blöbaum and Offen (2023) Ober-Blöbaum, S. and Offen, C. (2023).  Variational learning of Euler–Lagrange dynamics from data.  _Journal of Computational and Applied Mathematics_ , 421, 114780.  10.1016/j.cam.2022.114780.
  * Offen and Ober-Blöbaum (2022) Offen, C. and Ober-Blöbaum, S. (2022).  Symplectic integration of learned Hamiltonian systems.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 32(1), 013122.  10.1063/5.0065913.
  * Qin (2020) Qin, H. (2020).  Machine learning and serving of discrete field theories.  _Scientific Reports_ , 10(1), 1–15.
  * Saemundsson et al. (2020) Saemundsson, S., Terenin, A., Hofmann, K., and Deisenroth, M. (2020).  Variational integrator networks for physically structured embeddings.  In _International Conference on Artificial Intelligence and Statistics_ , 3078–3087. PMLR.
  * Santos et al. (2022) Santos, S., Ekal, M., and Ventura, R. (2022).  Symplectic momentum neural networks-using discrete variational mechanics as a prior in deep learning.  In _Learning for Dynamics and Control Conference_ , 584–595. PMLR.
  * Vermeeren (2017) Vermeeren, M. (2017).  Modified equations for variational integrators.  _Numerische Mathematik_ , 137(4), 1001–1037.


is a popular new peer-to-peer network protocol designed to
disseminate messages quickly and efficiently by allowing peers to
forward the full content of messages only to a dynamically selected
subset of their neighboring peers (mesh neighbors) while gossiping
about messages they have seen with the rest. Peers decide which of
their neighbors to graft or prune from their mesh locally and
periodically using a score for each neighbor. Scores are calculated
using a  that depends on mesh-specific parameters, weights
and counters relating to a peer's performance in the network. Since
a GossipSub network's performance ultimately depends on the
performance of its peers, an important question arises: Is the score
calculation mechanism effective in weeding out non-performing or
even intentionally misbehaving peers from meshes?  We answered this
question in the negative in our companion
paper <cit.> by reasoning about  using our
formal, official and executable ACL2s model. Based on our findings,
we synthesized and simulated attacks against  which were
confirmed by the developers of , , and , and publicly
disclosed in MITRE CVE-2022-47547. In this paper, we present a
detailed description of our model. We discuss design decisions,
security properties of , reasoning about the security properties
in context of our model, attack generation and lessons we learnt
when writing it.

§ INTRODUCTION

is a new peer-to-peer network protocol used by popular
applications like  <cit.> and
FileCoin <cit.>. Messages transmitted in a  network
are typically categorized into topics, which peers of the
network can subscribe to or unsubscribe from. A peer can be part of several
meshes corresponding to different topics. In contrast to flood
publishing where a peer forwards every message it receives to all
of its neighboring peers subscribed to the corresponding topic, uses lazy pull, wherein a peer forwards full messages only to
its mesh neighbors in the relevant topic. A peer can graft or
prune a mesh neighbor based on various heuristic security mechanisms
that ultimately rely on a locally computed score.  The score is
calculated periodically by each peer for each of its neighbors and is
never shared. The , which is used to calculate a neighboring
peer's score, depends on application-specific parameters and weights,
and takes into account the performance of the neighbor both generally
and on a given topic. Ideally, the score of misbehaving peers (,
peers that drop messages or forward invalid ones) is penalized, which
matters because negatively scored mesh neighbors get pruned.

The  developers specified their protocol in English
prose <cit.> and implemented it in
GoLang <cit.>. They relied on unit-tests and network
emulation <cit.> of pre-programmed
scenarios for testing to show that misbehaving peers in a  network
are eventually pruned. However, simple testing is not enough. Dijkstra
famously quoted: “Program testing can be a very effective way to show
the presence of bugs, but it is hopelessly inadequate for showing
their absence.”

In our companion paper <cit.>, we formalized the specification in the ACL2s (the ACL2
Sedan) <cit.> theorem prover. ACL2s extends
ACL2 <cit.> with an advanced data definition framework
(Defdata) <cit.>, the
cgen <cit.>
framework for automatic counterexample generation, a powerful
termination analysis based on calling context graphs <cit.> and
ordinals <cit.>, and a
property-based modelling/analysis framework, each of which helped
immensely in our formalization and verification effort. Our publicly
available model <cit.> is designed to be modular and pluggable
, it allows us to reason about parts of the model in isolation as
well as about applications running on top of the network. Officially,
does not come with any properties. We formalized and attempted to
prove security properties which (1) we thought should be reasonable
for a score-based protocol like  to satisfy and which (2) the developers agree with. One such property, which states that
continuously misbehaving peers are eventually pruned, turned out to be
invalid in the case of . We leveraged the cgen facility in ACL2s
to automatically discover vulnerable network states that invalidate
our property. We built attack gadgets using sequences of events which
can take a network from a reasonable starting state to one of the
discovered vulnerable states. We synthesized attacks which were
confirmed by the , , and  developers and publicly disclosed
in MITRE CVE-2022-47547.
At the time of writing, the  developers are actively working on a fix.
These results are explained in our companion
paper <cit.>. In this work, we focus on our
formalization of the  ACL2s model and its use for reasoning,
simulation and attack synthesis.

Paper Outline. Section <ref> describes the protocol and our ACL2s model
simultaneously. Section <ref> describes properties we
used to reason about . Based on the insights gleaned from
proving/disproving these properties, Section <ref>
describes how we synthesized attacks that can disrupt communication in
an   network. Section <ref> describes
some limitations of our model. Section <ref> presents related work
on mechanized theorem proving efforts in the field of distributed
systems. Section <ref> concludes.

§ GOSSIPSUB MODEL DESCRIPTION

In this section, we describe our ACL2s model, while also giving an
overview of how the  protocol works. Interested readers are
encouraged to walk through our ACL2s formalization whose presentation
mirrors this section <cit.>. Consider the mesh shown in
Figure <ref>.

A mesh of peers subscribing to the same topic.  This Figure was taken from <cit.>.

Full-message payloads are forwarded on the full-message peering edges,
which consume more bandwidth, while the metadata-only peering edges
are used only to advertise and request full-messages using
corresponding “IHAVE” and “IWANT” Remote
Procedure Calls (RPCs). These RPCs carry the metadata of the
full-message being advertised or requested, which is considerably
smaller than the message itself. In this way, network congestion is
kept in check, and full-messages are supposed to be eventually
disseminated to all peers who subscribe to the given topic.  As an
example, peers numbered 1 and 2 can send
full-messages to each other, since they are mesh neighbors on the
illustrated topic and share a full-message peering edge. However, peer
29 can receive a full-message from 1 only if it
requests one, by sending an IWANT message in response to a
corresponding IHAVE message received from peer 1.
Note, Figure <ref> only shows a single  mesh. However,
real world applications can have an arbitrary number of topics and
corresponding meshes, ,  supports up to 71
topics. Figure <ref> illustrates a community graph of
an actual  network from Li et. al. <cit.>.

A community graph of  nodes where a bigger node
implies greater degree. Similarly colored nodes have more edges
among them than to the rest of the graph.

In order to model and analyze  in ACL2s, we rely heavily on
Defdata to easily model network components, and on cgen for property-based testing.
Distributed systems are often described in terms
of state machines, namely automata that encode the discrete
behaviors of a peer in the system <cit.>.  In our model,
we describe the state-space of a  peer with a peer-state
type, and then implement the  state machine using a state
transition function. We use a Defdata map from peers to their
corresponding states to capture the state of an entire network.

(defdata group (map peer peer-state))

We use records whenever we need to store state because of the
convenience of using named fields to access the internals of the
state. While proving theorems referring to states, we noticed that we
routinely had to prove helper lemmas about the types of the contents
of those states. This motivated us to rewrite records in the ACL2s
books to enable such helper lemmas automatically, leading
to cleaner code.

The local state of a peer is modeled as a record using the following definition:

(defdata peer-state
(record (nts        . nbr-topic-state)
(mst        . msgs-state)
(nbr-tctrs  . pt-tctrs-map)
(nbr-gctrs  . p-gctrs-map)
(nbr-scores . peer-rational-map)))

where (1) nts is a record that stores information about the
peer's neighbors' subscriptions, the peer's mesh neighbors, and the
peer's fanout (which we describe later); (2) mst is a record
that stores the peer's messages and related state; (3)
nbr-tctrs and nbr-gctrs are total maps that store
counters used for computing a neighbor's topic-specific and general
scores (respectively); and finally (4) nbr-scores is a total
map from peers to their scores. Note that peer and
topic are ACL2 symbols, while tctrs is a record that
keeps track of a peer's behaviors in each topic:

(defdata tctrs
(record (invalidMessageDeliveries . non-neg-rational)
(meshMessageDeliveries    . non-neg-rational)
(meshTime                 . non-neg-rational)
(firstMessageDeliveries   . non-neg-rational)
(meshFailurePenalty       . non-neg-rational)))

Similarly, gctrs keeps track of a peer's general behaviors, not
pertaining to any single topic:

(defdata gctrs
(record (apco . rational)  ;; application provided score
(ipco . non-neg-rational) ;; ip-colocation factor
(bhvo . non-neg-rational))) ;; track misbehavior

Notice that each of the counters is a rational, not a natural. This is
because counters are supposed to fractionally decay at regular
intervals. The rate of decay is dependent on the application running
on top of .  We define maps for each of these counters:

(defdata pt (cons peer topic))
(defdata pt-tctrs-map (map pt tctrs))
(defdata p-gctrs-map (map peer gctrs))
(defdata peer-rational-map (map peer rational))

For each map, we define a lookup function with default values.
automatically proves termination and input-output contracts,
which suffices to show that the maps are total.
As an example, the following is the lookup function for scores:

(definecd lookup-score (p :peer prmap :peer-rational-map) :rational
(let ((x (mget p prmap)))
(match x
(nil 0)  ;; default value
(    x))))

We now explore each component of the peer-state.

Neighbor topics state. A key objective of the protocol is to reduce network congestion, while forwarding data as
quickly as possible. To achieve this, a peer forwards full messages
only to a subset of its neighboring peers. A  peer keeps
track of the topics its neighbors subscribe to, using the
nbr-topicsubs field in nbr-topic-state. Using this
information, it is able to build up a picture of the topics around it
and which peers are subscribed to each topic. The peers to which it
forwards full messages in topics it does not itself subscribe to constitute
a list called a fanout, and is stored in the topic-fanout field,
which is a map from topics to lists of peers (lop). The peer
forwards full messages to other peers with whom it shares mesh membership.
Mesh memberships are stored in the topic-lop-map
field topic-mesh. Finally, the peer stores the last time it
published a message to its fanout. This is used to expire a peer's
fanout, if it has been too long since the peer last published.

(defdata peer symbol)
(defdata lop (listof peer))
(defdata topic-lop-map (map topic lop))
(defdata topic-nnr-map (map topic non-neg-rational))

(defdata nbr-topic-state
(record (nbr-topicsubs . topic-lop-map)
(topic-fanout  . topic-lop-map)
(last-pub      . topic-nnr-map)
(topic-mesh    . topic-lop-map)))

Messages state. msgs-state is a record type used to
store information about messages received, requested, seen, or
forwarded by a peer. First we define the types pid-type,
which is just an alias for the type symbol, and the record
type payload-type.

(defdata-alias pid-type symbol)
(defdata payload-type (record (content . symbol)
(pid     . pid-type)
(top     . topic)
(origin  . peer)))

payload-type is a record used to represent full messages, carrying
the message content, payload id, the topic of this message and the
peer who originated it. pid-type represents payload ids,
a hash of a message payload which can identify the full
message content. The msgs-state is defined as follows:

(defdata msg-peer (v (cons payload-type peer)
(cons pid-type peer)))
(defdata msgpeer-rat (map msg-peer rational))
(defdata msgs-waiting-for (map pid-type peer))
(defdata mcache (alistof payload-type peer))

(defdata msgs-state
(record (recently-seen   . msgpeer-rat)
(pld-cache       . mcache)
(hwindows        . lon)
(waitingfor      . msgs-waiting-for)
(served          . msgpeer-rat)
(ihaves-received . nat)
(ihaves-sent     . nat)))

msgs-state stores (1) recently-seen, a map from
either a full-message or a message hash to the time since
receipt; (2) pld-cache, an association
list of full messages and their senders; (3) hwindows, or
history windows, a list of naturals where each is the
number of messages received in an interval; (4) waitingfor, a
map from message ids of messages that haven't been received yet, to
peers one has sent corresponding IWANT requests to; and
(5) served, a map from
either a full-message or a message hash to a rational, denoting the
count of the number of times this message was served. The served map
helps to
detect peers sending too many IWANT messages.
Finally, msgs-state contains (6)
ihaves-received and ihaves-sent, which store the number of
IHAVE messages received and sent, respectively.

Events. Our model includes events that can occur in a
network. Events include peers sending or receiving (1) control
messages like GRAFT or PRUNE for mesh control,
SUB, UNSUB, JOIN and LEAVE for
topics, or CONNECT1 or CONNECT2 messages to edit
neighbors; (2) PAYLOAD for carrying message payload, or
IHAVE for advertising or IWANT for requesting
message payloads; and (3) HBM for heart-beat events occurring
at each peer at regular intervals. A list of events will have type
loev. Events are defined as follows:

(defdata verb (enum '(SND RCV)))
(defdata rpc (v (list 'CONNECT1 lot)
(list 'CONNECT2 lot)
(list 'PRUNE topic)
(list 'GRAFT topic)
(list 'SUB topic)
(list 'UNSUB topic)))
(defdata data (v (list 'IHAVE lopid)
(list 'IWANT lopid)
(list 'PAYLOAD payload-type)))
(defdata mssg (v rpc data))
(defdata evnt (v (cons peer (cons verb (cons peer mssg)))
(list peer 'JOIN topic)
(list peer 'LEAVE topic)
(list peer verb peer 'CONNECT1 lot)
(list peer verb peer 'CONNECT2 lot)
(list peer 'HBM pos-rational)
(list peer 'APP payload-type)))
(defdata hbm-evnt (list peer 'HBM pos-rational))
(defdata-subtype hbm-evnt evnt)

(defdata loev (listof evnt))

Heart-beat events occur at regular intervals at each peer in a network. Several maintenance activities are performed during each such
event. Scores are updated for neighboring peers, which are then used
to update mesh memberships. Counters are multiplied by some
application-specific decay factors.  Neighboring peers that are not
part of any mesh are sent GRAFT messages at regular
intervals, provided that their addition could improve the average
score of peers in the corresponding mesh. Several of these actions
depend on application-specific weights (used for calculating scores)
and parameters.

Parameters and Scoring. We store the weights and parameters
relevant for scoring in a record called a twp.
This record totally captures the application-specific configuration
of a  instance, , we can uniquely specify the configuration used by in a twp.  Thus, in order to simulate an application running on top of
with our model, all we need to know is the application-specific twp.
This is what we mean when we say our model is “pluggable”.

(defdata weights
(record (w1  . non-neg-rational)
(w2  . non-neg-rational)
(w3  . non-pos-rational)
(w3b . non-pos-rational)
(w4  . neg-rational)
(w5  . non-neg-rational)
(w6  . neg-rational)
(w7  . neg-rational)))

(defdata params
(record (activationWindow               . nat)
(meshTimeQuantum                . pos)
(p2cap                          . nat)
(timeQuantaInMeshCap            . nat)
(meshMessageDeliveriesCap       . pos-rational)
(meshMessageDeliveriesThreshold . pos-rational)
(topiccap                       . rational)
(grayListThreshold              . rational)
(d                              . nat)
(dlow                           . nat)
(dhigh                          . nat)
(dlazy                          . nat)
(hbmInterval                    . pos-rational)
(fanoutTTL                      . pos-rational)
(mcacheLen                      . pos)
(mcacheGsp                      . non-neg-rational)
(seenTTL                        . non-neg-rational)
(opportunisticGraftThreshold    . non-neg-rational)
(topicWeight                    . non-neg-rational)
(meshMessageDeliveriesDecay     . frac)
(firstMessageDeliveriesDecay    . frac)
(behaviourPenaltyDecay          . frac)
(meshFailurePenaltyDecay        . frac)
(invalidMessageDeliveriesDecay  . frac)
(decayToZero                    . frac)
(decayInterval                  . pos-rational)))

(defdata wp (cons weights params))
(defdata twp (map topic wp))

Note that  required weights w3 and w3b to be
negative. However, the use of Defdata allowed us to automatically find
that  used an invalid twp because it set w3 and
w3b to zero. We discussed this with the  developers who
then agreed to allow zero values.  Given T, a set of topics
which the neighbor subscribes to; tctrs, the neighbor's topic
specific counters; gctrs, the neighbor's global counters; and
a twp containing entries for each topic our neighbor
subscribes to, the  calculates a neighbor's score as shown

\[
\begin{aligned}
\mathit{score}(q)
& = \min\big(\text{TC},\sum_{\tau \in T} \text{tw}^\tau \times
\sum_{i \in \{1,2,3,3b,4\}} w_i^\tau P_i^\tau \big) +
w_5 P_5 +
w_6 P_6 +
w_7 P_7
\end{aligned}
\]

\begin{align*}
\codify{(weights . params)} &= \codify{(mget $\tau$ twp)}\\
w_i^\tau &= \codify{(weights-w\emph{i} weights)}\\
P_1^\tau &= \codify{(calcP1}\codify{(tctrs-meshTime tctrs) (params-meshTimeQuantum params)}\\
&\quad\ \codify{(params-timeQuantaInMeshCap params))} \\
P_2^\tau &= \codify{(calcP2}\codify{(tctrs-firstMessageDeliveries tctrs) (params-p2cap params))}\\
P_3^\tau &= \codify{(calcP3}\codify{(tctrs-meshTime tctrs) (params-activationWindow params)} \\
&\quad\ \codify{(tctrs-meshMessageDeliveries tctrs)}\\
&\quad\ \codify{(params-meshMessageDeliveriesCap params)}\\
&\quad\ \codify{(params-meshMessageDeliveriesThreshold params))}\\
P_{3b}^\tau &= \codify{(calcP3b}\codify{(tctrs-meshTime tctrs) (params-activationWindow params)}\\
&\quad\ \codify{(tctrs-meshFailurePenalty tctrs)}\\
&\quad\ \codify{(tctrs-meshMessageDeliveries tctrs)}\\
&\quad\ \codify{(params-meshMessageDeliveriesCap
&\quad\ \codify{(params-meshMessageDeliveriesThreshold params))}\\
P_4^\tau &= \codify{(calcP4}\codify{(tctrs-invalidMessageDeliveries
tctrs))} \\
P_5 &= \codify{(gctrs-apco gctrs)} \\
P_6 &= \codify{(gctrs-ipco gctrs)} \\
P_7 &= \codify{(calcP7 (gctrs-bhvo gctrs))}\\
\text{tw}^\tau &= \codify{(params-topicweight params)}\\
\text{TC} &= \codify{(params-topiccap (cddar twp))}
\end{align*}

TC does not depend on any topic, but since it is stored in a
twp which is indexed by topic, its value is
replicated in each of the corresponding params. So, it is
fine to extract its value from the first entry of a twp. Note
that for score calculations, we require a non-empty twp.

Each of the calcPi functions where
$\codify{i} \in \{\codify{1,2,3,3b,4,7}\}$ is used to calculate
contributions to the score by one or more of counter values from
tctrs. calcP1 calculates the contribution to a neighbor's
score based on the time spent in common meshes. calcP2 awards score
for being one of the first few to forward a message. calcP3
calculates penalties due to the mesh message transmission rate being
below a given threshold of
(params-meshMessageDeliveriesThreshold params). Whenever a
peer is pruned, its corresponding tctrs counter
meshFailurePenalty is augmented by the mesh message
transmission rate deficit. This counter is not cleared even after the
peer has been pruned. calcP3b scores mesh message delivery failures
based on the value of this counter. Hence, $P^\tau_{3b}$ is a “sticky”
value which is supposed to discourage a peer that was pruned because
of under-delivery from quickly getting re-grafted in a mesh. calcP4
calculates the penalty on score due to sending invalid
messages. calcP7 calculates penalties due to several kinds of
misbehaviors described by the  specification. An example of such
misbehaviors includes spamming with too many IHAVE messages
which are either bogus and/or not following up to the corresponding
IWANT requests.

The Transition Function. We define a transition function
run-network, which, given an initial Group state and
a list of evnt, produces a trace of type egl, which
is an alist of evnt and group. Hence, after running
a simulation, we have access to the state of the Group after
each evnt was processed. run-network depends on the
transition function for the peer-state (transition),
which depends on the transition functions for nbr-topic-state
(update-nbr-topic-state) and for msgs-state

(update-msgs-state). For brevity, we mention only the
signatures of each of these functions below. Notice that each of these
signatures represents neatly and concisely the types of the formal
arguments and the function return type, which is very useful in a
large code base.

(defdata egl (alistof evnt group)) ;; simulation trace

(definecd run-network (gr :group evnts :loev i :nat r :twp s :nat) :egl

(defdata peer-state-ret
(record (pst . peer-state)
(evs . loev)))

(definecd transition
(self :peer pstate :peer-state evnt :evnt r :twp s :nat) :peer-state-ret

(defdata msgs-state-ret
(record (mst . msgs-state)
(evs . loev)
(tcm . pt-tctrs-map)
(gcm . p-gctrs-map)))

(definecd update-msgs-state (mst :msgs-state evnt :evnt pcm :pt-tctrs-map
gcm :p-gctrs-map r :twp) :msgs-state-ret

(defdata nbr-topic-state-ret
(record (nts . nbr-topic-state)
(evs . loev)
(tcm . pt-tctrs-map)
(gcm . p-gctrs-map)
(sc  . peer-rational-map)))

(definecd update-nbr-topic-state (nts :nbr-topic-state
nbr-scores :peer-rational-map
tcm :pt-tctrs-map gcm :p-gctrs-map
evnt :evnt r :twp s :nat) :nbr-topic-state-ret

A  peer can select a random subset of its fanout and promote them
as mesh members. It can also selects a random subset of its neighbors
to advertise with “IHAVE” messages. Such non-determinism is handled
by sending a random seed s as a formal parameter to
run-network, which is then propagated to the other transition
functions it depends on. Observe that a single event like full message
forwarding can trigger several more forwards, causing a cascade of
events. Such events are represented by the evs field in the
return types of update-msgs-state and
update-nbr-topic-state. In order to limit the total number of
events processed, we send a natural number i as a formal
parameter to the run-network function.

When proving contract theorems for the transition functions, we needed
to prove the types of terms returned by utility functions, like
shuffle, or ACL2 functions like
set-difference-equal. For this, we used polymorphism and
automated type-based reasoning provided by Defdata, as shown below:

(sig set-difference-equal ((listof :a) (listof :a)) => (listof :a))
(sig shuffle ((listof :a) nat) => (listof :a))

We made heavy use of higher-order macros written by
Manolios <cit.> to improve the readability of our
code. For example, create-map* is a list functor. Given an
admitted function name or a lambda expression f of type
a $\rightarrow$ b, create-map* defines a function
map*-*f of type (listof a) $\rightarrow$ (listof
b). In the following code snippet, we show theorems proving that it
obeys the functor laws, and give an example of its usage. map* is a
syntactic sugar that maps function f onto a list without referring to the
generated function name map*-*f.

(definec id (x :all) :all
;; Proof that create -map* is a list functor
;; 1) functor mapping preserves the identity function
(property functor-id (xs :tl)
(== (map* id xs) ;; id is an identity function for lists
(id xs)))

;; 2) functor mapping preserves function composition. f and g are declared
;; using defstub. gof is defined as a composition of g and f
(property functor-comp (xs :tl)
(== (map* gof xs)
(map* g (map* f xs))))

;; function to create a list of SND GRAFT events from peer p to a list of peers
(create-map* (lambda (tp p) `(,p SND ,(cdr tp) GRAFT ,(car tp)))
(:name mk-grafts)
(:fixed-vars ((peerp p))))

(check= (map* mk-grafts '((FM . A) (DS . B)) 'P)
'((P SND A GRAFT FM) (P SND B GRAFT DS)))

Given an admitted function name or a lambda expression f of
type a $\times$ b $\rightarrow$ b, the higher order function
create-reduce* defines a function reduce*-*f which
accepts a list of elements of type a, an initial accumulator
value of type b, and returns a reduction of the list using
f, from left to right.

;; function to extract all the subscribers (neighboring peers) from a topic-lop-map
(create-reduce* (lambda (tp-ps tmp) (app tmp (cdr tp-ps)))
(:name subscribers))

(check= (reduce* subscribers '()
'((T1 P1 P2 P3)
(T2 P4 P5 P1)))
'(P4 P5 P1 P1 P2 P3))

§ REASONING ABOUT THE SCORING FUNCTION

Based on the observation that honest peers can be distinguished from
malicious ones based on their observable behaviors (using local
counters and scores), and thus, the overall network can be made more
secure and performant if every honest peer promotes their
well-behaving neighbors and demotes poorly-behaved ones, we came up
with the following informal fundamental property.

Fundamental Property of  Defense Mechanisms. Peers
who behave poorly will be demoted by their neighbors.  Peers who
behave better-than-average will be promoted by their neighbors.
Promotion/demotion is entirely based on local peer behavior.

An example network.

Before reasoning about the formalization of this fundamental property
later in the paper, we discuss the importance of this fundamental
property. Consider the simple example shown in
Figure <ref>, where peers A and B subscribe to, and are
mesh neighbors in both Red and Blue topics. It might be possible for B
to observe that A behaves perfectly well in the Blue topic while
simultaneously misbehaving in the Red topic. This is not good for B
because it depends on A for all of its messages in the Red topic. Note
that this is a very simplified example. In an actual network, B could
have several other neighbors subscribed to the Red topic. But as we will
show later, it is equally trivial to have a scenario where all of B's
neighbors isolate it from communications in the Red topic. In this
example, we want B to prune A from its Red mesh in hopes of finding a
better mesh neighbor later on. Reasoning about this fundamental
property directly would be difficult due to the massive search-space
of possible attack vectors. Hence, we focus on the following liveness
property capturing the essence of the fundamental property:

If a peer's score relating to its performance in any topic is
continuously non-positive, then the peer's overall score should
eventually be non-positive:

\[\forall q,\tau::\langle
\glb(\mathit{score(q)}\textit{ for topic }\tau \leq 0) \Rightarrow
\fin(\mathit{score(q)} \leq 0) \rangle \]

where score$(q)$ for topic $t$ is defined below.

\[
\mbox{tw}^\tau \times \sum_{i \in \{1,2,3,3b,4\}} w_i^\tau P_i^\tau
\]

Notice that Property <ref> is temporal. We write the non-temporal
version of this property in context of  (using twp) in ACL2s as shown below, and disprove the temporal version
using an induction argument later in the paper.

(property (ptc :pt-tctrs-map pcm :p-gctrs-map p :peer top :topic)
:hyps  (^ (member-equal `(,p . ,top) (acl2::alist-keys ptc))
(> (lookup-score p (calc-nbr-scores-map ptc pcm *eth-twp*)) 0))
(> (calcScoreTopic (lookup-tctrs p top ptc) (mget top *eth-twp*)) 0))

The following is one of the counter-examples to the above property,
generated by cgen in :

((top 'agg)
(p 'p4)
(pcm '((p3449 (:0tag . gctrs) (:apco . 0) (:bhvo . 0) (:ipco . 0))
(p3450 (:0tag . gctrs) (:apco . 0) (:bhvo . 0) (:ipco . 0))
(p3451 (:0tag . gctrs) (:apco . 0) (:bhvo . 0) (:ipco . 0))))
(ptc '(((p4 . agg)
(:0tag . tctrs)
(:firstmessagedeliveries . 0)
(:invalidmessagedeliveries . 0)
(:meshfailurepenalty . 0)
(:meshmessagedeliveries . 1)
(:meshtime . 42))
((p4 . blocks)
(:0tag . tctrs)
(:firstmessagedeliveries . 324)
(:invalidmessagedeliveries . 0)
(:meshfailurepenalty . 0)
(:meshmessagedeliveries . 330)
(:meshtime . 377))
((p4 . sub1)
(:0tag . tctrs)
(:firstmessagedeliveries . 371)
(:invalidmessagedeliveries . 0)
(:meshfailurepenalty . 0)
(:meshmessagedeliveries . 377)
(:meshtime . 324))
((p4 . sub2)
(:0tag . tctrs)
(:firstmessagedeliveries . 318)
(:invalidmessagedeliveries . 0)
(:meshfailurepenalty . 0)
(:meshmessagedeliveries . 324)
(:meshtime . 371))
... ))

For brevity, we omit entries for peer-topic key values for peers other
than p4. In property-based testing, the free variables of a property under test
are assigned values using a synergistic combination of theorem proving
and random assignments computed using type-based enumerators
(generators) in an effort to discover counterexamples to the property.
Observe that Property <ref> depends on ptc and
pcm which do not have trivial types. These are maps
containing records which themselves consist of several numerical
values, which makes the search space of possible counter-examples
immensely large. In order to make it easier for cgen to find
counter-examples, we wrote custom enumerators to enumerate restricted
values for topic, tctrs, gctrs,
pt-tctrs-map and p-gctrs-map. Specifically, we limit
the penalties on the score due to some counters, such that the
negative contributions due to misbehavior are comparable to the
positive contributions due to good behavior. Below we define custom
enumerators for topic and tctrs.

(definec topics () :tl
;; valid topics used in Ethereum
'(AGG BLOCKS SUB1 SUB2 SUB3))

(definec nth-topic-custom (n :nat) :symbol
(nth (mod n (len (topics))) (topics)))
(defdata lows (range integer (0 <= _ <= 1))) ;; high values
(defdata-subtype lows nat)
(defdata highs (range integer (300 < _ <= 400))) ;; low values
(defdata-subtype highs nat)

(defun nth-bad-counters-custom (n)
;; setting invalidMessageDeliveries and meshFailurePenalty to 0 due to high penalty
(defun nth-good-counters-custom (n)
(tctrs 0 (nth-highs (+ n 2)) (nth-highs (+ n 3)) (nth-highs (+ n 4)) 0))

(defun nth-counters-custom (n) ;; custom enumerator for tctrs
(if (== 0 (mod n 4))
(nth-bad-counters-custom n)
(nth-good-counters-custom n)))

Besides Property <ref>, we formalize three safety properties
for GossipSub, stated below. Together, these four are the most
general properties of the score function, which must hold in order
for the fundamental property to hold.

Increasing bad-performance counters (which are multiplied with
negative weights) should decrease the overall score.

Increasing good-performance counters (which are multiplied with
positive weights) will not decrease the score for a mesh peer that has
been in the mesh for a sufficiently long time.

If two peers subscribe to the same topics, and achieve identical
per-topic counters, and identical global counters, then they achieve
identical scores.

We are able to find counterexamples to Property <ref> in much the
same way as <ref>, however, in this work we focus on the
counterexamples to Property <ref>, which are more interesting
in terms of attack generation.  We manually prove that
Property <ref> holds over all configurations (available with the
paper artifacts).  The proof that Property <ref> holds over all
configurations follows directly from referential transparency.

We also prove a limit on the maximum score achievable in a topic, as
shown below.
(property max-topic-score (tctrs :tctrs weights :weights params :params)
(<= (calcScoreTopic tctrs (cons weights params))
(* (params-topicweight params)
(+ (* (mget :w1 weights) (params-timeQuantaInMeshCap params))
(* (mget :w2 weights) (params-p2cap params))))))

§ ATTACK GENERATION

The counter-example <ref> which we obtained in the previous
section is a specification for an unsafe state that does not satisfy
Property <ref> for an  network. However, we are also
interested in characterizing and generating attacks against an  network for three main reasons: (1) to show that an unsafe state
is reachable from a reasonable start state, (2)
invalidation  of our temporal Property <ref> using the
trace generated from the attack, and finally (3) demonstration of
scalability of our attack to large networks, typically of the
size and shape used by real world applications.

Counter-example <ref> suggests that an unsafe state is one where
a neighboring peer throttles communication in a particular topic
while maintaining an overall positive score, hence, avoiding getting
pruned. Using this insight, we design attack gadgets that can
perpetrate such attacks locally. We define an attack gadget as a tuple
$\langle A, V, S \rangle$, where $A, V$ are peers ($A$ is the attacker
and $V$ is the victim), $S$ is a set of subnet topics (the attacked
topics), and $A, V$ are mesh neighbors over a set of topics that is a
superset of $S$. For each $i \in \mathds{N}$, we define i to be
the set of attack gadgets where $|S|=i$. Figure <ref>
illustrates an example attack gadget in 2 where peers $A$ and $V$
are neighbors in four meshes corresponding to topics: Red, Yellow,
Blue and Green, out of which $A$ is attacking $V$ in
$S=\{\mbox{Red}, \mbox{Blue}\}$.

An example 2 attack gadget $\langle A, V, \{\mbox{Red}, \mbox{Blue}\} \rangle$.

We generate a sequence of events consisting of message transmission
events from $A$ to $V$ (referred to as a and v in
code) as well as heart-beat events at $V$ ($V$ updates the scores of
its peers during heart–beat events). These events are designed to
either restrict or completely block communication to $V$ in the
attacked topics while maintaining normal communication rate in all the
other topics, as shown in the following code snippet.

(definecd emit-evnts (a v :peer ts ats :lot n m e :nat) :loev
;; mesh message deliveries in attacked topics
(app (emit-meshmsgdeliveries-peer-topics a v ats m)
;; mesh message deliveries in other topics
(emit-meshmsgdeliveries-peer-topics a v (set-difference-equal ts ats) n)
;; heart-beat events at the victim node
`((,p2 HBM ,e))))

The expression (emit-meshmsgdeliveries-peer-topics a v ts ats
n m e) generates a list of events $E$ sending m mesh
messages in the attacked topics and n mesh messages in the
other topics from peer a to peer v per heart-beat at
v which happens every e seconds. m is
generally set to 0 or 1 in order to block or
throttle communication in the attacked topic meshes. We choose a start
state of the network based on actual full-network topologies of the
testnets Ropsten (shown in Figure <ref>), Goerli and
Rinkeby, as measured by Li
et. al <cit.>. Table <ref> characterizes each
of these topologies.

min     max     avg

Ropsten      588       1      418       25.49       5
Goerli       1355      1      712       28.26      5
Rinkeby      446       1      191       68.96      6

Network Characteristics

We create our start state as a Group, using the topologies
provided. In this Group, we initialize our attack gadget and
simulate its run over the generated sequence of events $E$ using the
run-network function. We use function
scorePropViolation described below, to
detect violations of Property <ref> in peer states occurring
in the output trace.

;; ats is the list of topics being attacked
(definec scorePropViolation (ps :peer-state p :peer ats :lot twpm :twp) :boolean
(match ats
(() (> (lookup-score p (calc-nbr-scores-map (peer-state-nbr-tctrs ps)
(peer-state-nbr-gctrs ps) twpm))
((top . rst) (^ (< (calcScoreTopic
(lookup-tctrs p top (peer-state-nbr-tctrs ps))
(mget top twpm))
(scorePropViolation ps p rst twpm)))))

We observe that after the first heart-beat event at the victim peer,
the state of the network at each subsequent heart-beat is
identical. Since Property <ref> is being violated at each of
these events, we make an inductive claim that it will forever be
violated, thus giving a counter-example to our liveness
Property <ref>. We wrote optimized versions of the
run-network function to generate traces of property
violations as a list of booleans, instead of generating the full trace
of type egl. We ran our simulations for 100,000 events to
ensure that we ended up in identical states, taking about a
minute on an M1 Macbook Air.

Finally, we scaled our attacks to build eclipse and
network partition attacks using a combination of attack
gadgets. <Ref> give
an intuition of how to combine our attack gadgets to carry out these
attacks. Our companion paper discusses the specifics of these attacks
in more detail.

An eclipse attack using $AG_2$ gadgets

A partition attack using $AG_2$ gadgets

§ LIMITATIONS

We now discuss limitations of our model. The most
crucial aspect of property-based testing is counter-example generation for invalid
properties. As explained previously, our properties depend on
complex types, for which we had to write custom
enumerators. Coming up with enumerators that had a high
probability of satisfying the hypotheses of our properties required considerable
analysis of  so as to
restrict certain tctrs values from skewing the
scores too much. Testing properties for new applications will likewise require
writing new custom enumerators. One might need to write a new
event generator as well, possibly generating events of shapes different
from the ones we showed.

§ RELATED WORK

The Protocol Labs ResNetLab and software audit firm Least Authority
tested  against a list of specific pre-programmed attack
scenarios <cit.> designed to degrade overall network
performance using a network emulator called
Testground <cit.>. Due to their use of simplified
configurations with only one topic (and because simple testing is not
enough to find bugs) they found that all the attacks failed against
, and the  made  more resilient to malicious nodes
attacks than the other tested protocols <cit.>.
Separate from simulation testing, Least Authority also audited the
implementation and provided recommendations for
improvement <cit.>. Though our work contributes the first
formalization of , there has been considerable previous work on
utilizing formal methods to
reason about distributed systems. We survey such works below.

Model Checking based approaches.  Lamport's modeling language
TLA+ <cit.> and the corresponding TLC
model checker <cit.> have been used to analyze properties
of distributed systems including Disk
Paxos <cit.>,
MongoRaftReconfig <cit.>, Byzantine
Paxos <cit.>,
Spire <cit.>, etc. McMillan and Zuck
applied specification-based testing to the Quic protocol, and
found vulnerabilities <cit.>.  Wu et. al. formally
modeled the Bluetooth stack using ProVerif, a model checker,
and found five known vulnerabilities and two new
ones <cit.>.  Chothia et. al. demonstrated the use
of ProVerif to verify distance-bounding protocols, , those
used by MasterCard and NXP <cit.>.  Separately,
Chothia modeled the MUTE anonymous file-sharing system using
the $\pi$-calculus, and proved the system insecure (discovering a novel
attack) <cit.>.  Cremers et. al. modeled all
handshake modes of TLS 1.3 using Tamarin, another model
checking tool, and discovered an unexpected
behavior <cit.>. An issue with using model checking
tools like ProVerif or Tamarin to verify a protocol like  is the
immense size of the state space needed to be checked, making them
infeasible for our use.

Refinement-based proof formalization. The theory of
refinement has proved to be useful for enabling the mechanical
verification of distributed systems' properties. Manolios' work on
refinement <cit.> has been
previously used for mechanical verification of pipelined
processors <cit.>. Manolios et. al. combined
theorem proving (using refinement maps) with model checking to verify
the alternating bit
protocol. <cit.>. IronFleet <cit.>
refines TLA style state-machine specification of a
Paxos-based library and a sharded key-value store to low
level implementation in Dafny (a SMT based program verifier) for
verification using Hoare-logic. Woo et. al. <cit.>
formally verified 90 properties of the RAFT protocol using
Verdi <cit.>, a tool they built in the
Coq proof assistant. Though they did not build an executable
model, their framework can be used to extract an executable protocol
implementation in OCaml. Verdi provides verified system
transformers used to refine a system in an ideal fault model to a
more realistic fault model, without any proof overhead on part of the
user. We believe that looking at  through the lens of refinement
will be interesting because, not only will it allow us to explain why
it failed our properties, but also guide us towards improving it.

Inductive-invariant based proof formalization.
Padon et. al. <cit.> proved the correctness of a simple
model of Paxos described in Effectively Propositional Logic (a
decidable fragment of First Order Logic) using
IVY <cit.>, a SMT-based safety verification
tool.  IVY can be used for verifying inductive invariants
about global states of a distributed protocol. Both the modeling and
the specification languages of IVY are restricted to a
decidable fragment of First Order Logic to ensure that all
verification conditions can be checked algorithmically. Hippel
et. al. <cit.> also used IVY to formally describe and
reason about Karn's Algorithm, a mechanism used to study rount trip
times of message transmissions. However, since IVY lacks a theory of
rationals, modelling the scoring function of  would not have been
possible using this tool.

Full stack verification.
Certain high-assurance distributed systems might require the whole
stack to be formally-verified. Such applications could, for instance,
be implemented on top of sel4: a high-performance operating
system microkernel that was formally verified against an abstract
specification using higher-order logic <cit.>. Another
example is the fully verified CLI stack <cit.>, a system
comprising of an operating system with some applications running in
it, operational semantics for two high level languages, a stack based
assembly language and the instruction set architecture (ISA), all the
way down to the register transfer level (RTL) design for a
microprocessor. The full stack was verified in Nqthm <cit.>.

§ CONCLUSION AND FUTURE WORK

In this paper, we described the  protocol, as well as our official
formalization based on its prose specification using the  theorem
prover. We explained our state models, transition functions as well as
design decisions. We showed our security property for  and how we
were able to find counter-examples against it. Finally, we described
several kinds of attacks we synthesized based on our attack gadgets,
using the couter-examples as specifications.

In the future, we would like to characterize an ideal variant of as a refinement of simpler protocols so as to prove safety properties,
as well as to contrast the ideal variant with our current model in
order to better explain why it is susceptible to attacks from
misbehaving peers. We would also like to support reasoning for the
application layer on top of our network model layer, since interesting
bugs can be found at the interface of these two layers.