A Riemannian Framework for Learning Reduced-order Lagrangian Dynamics

## Abstract

Abstract By incorporating physical consistency as inductive bias, deep neural networks display increased generalization capabilities and data efficiency in learning nonlinear dynamic models. However, the complexity of these models generally increases with the system dimensionality, requiring larger datasets, more complex deep networks, and significant computational effort.
We propose a novel geometric network architecture to learn physically-consistent reduced-order dynamic parameters that accurately describe the original high-dimensional system behavior.
This is achieved by building on recent advances in model-order reduction and by adopting a Riemannian perspective to jointly learn a structure-preserving latent space and the associated low-dimensional dynamics.
Our approach enables accurate long-term predictions of the high-dimensional dynamics of rigid and deformable systems with increased data efficiency by inferring interpretable and physically plausible reduced Lagrangian models.

## 1 Introduction

Deep learning models recently emerged as powerful tools for learning the continuous-time dynamics of physical systems. In contrast to classical system identification methods which rely on analytical derivations, black-box approaches learn to predict the behavior of nonlinear systems using large amounts of trajectory data.
However, collecting such data is prohibitively expensive, and the learned models often predict trajectories that violate the laws of physics, e.g., by not conserving energy.

Gray-box approaches tackle such limitations by incorporating physics priors as inductive biases into deep learning architectures. These methods have been broadly investigated for different systems, such as Lagrangian and Hamiltonian mechanics, and function approximators, including neural networks (Lutter et al., ; Greydanus et al., ; Cranmer et al., ; Lutter & Peters, ) and Gaussian processes (Tanaka et al., ; Evangelisti & Hirche, ).
Several physically-consistent models were extended to capture non-conservative forces and contact mechanics (Hochlehnert et al., ; Zhong et al., ), and to ensure longer-term stability of the identified dynamics by including differential solvers (Chen et al., ) into the training loop (Finzi et al., ; Zhong & Leonard, ). The model performances were also improved by representing the dynamics using coordinates adapted to the physical system at hand (Finzi et al., ; Celledoni et al., ; Duong & Atanasov, ). Altogether, gray-box approaches are physically consistent, versatile, more data-efficient, and display increased generalization capabilities compared to black-box methods (Lutter & Peters, ; Greydanus et al., ).
Still, gray-box approaches have, so far, only managed to learn the dynamics of low-dimensional physical systems, leaving the question of modeling high-dimensional dynamical systems, such as fluid flows, continuum mechanics, and robots, unanswered.

Predicting trajectories of high-dimensional systems is notoriously difficult due to the computational cost of solving high-dimensional and highly nonlinear differential equations. In this context, model order reduction (MOR) techniques find a computationally efficient yet descriptive low-dimensional surrogate system — a reduced-order model (ROM) — of a given high-dimensional dynamical system or full-order model (FOM) with known dynamics (Schilders et al., ). MOR is typically achieved by projecting the FOM onto a lower-dimensional space via linear (Thieffry et al., ; Farhat et al., ; Carlberg et al., ) or nonlinear (Sharma et al., ; Barnett & Farhat, ) mappings.
In particular, Autoencoders were shown to be well suited to extract descriptive nonlinear reduced representations from data (Lee & Carlberg, ; Buchfink et al., ; Otto et al., ).
Recent approaches investigated the preservation of physics-induced geometric structures, i.e. Lagrangian or Hamiltonian structures, of FOMs within projection-based ROMs (Carlberg et al., ; Lee & Carlberg, ; Otto et al., ; Hesthaven et al., ). This structure awareness enables the derivation of shared theoretical properties such as energy conservation and stability preservation (Buchfink et al., ), as well as the design of low-dimensional control strategies on the ROM (Lepri et al., ).
Buchfink et al. () recently unified various structure-preserving MOR techniques by adopting a differential-geometric perspective on the problem. In this framework, ROMs are defined on embedded submanifolds of the high-dimensional Riemannian or symplectic manifold associated with a given Lagrangian or Hamiltonian system.

Arguably the complexity of learning the dynamics of physical systems scales with their dimensionality, as the complexity of dynamic behaviors generally increases and requires more complex network architectures and more training data.
*In this paper*, we propose a novel physics-inspired geometric model that learns the continuous-time dynamics of high-dimensional systems.
Our *first contribution* is to adopt a differential geometry perspective and jointly learn *(1)* a structure-preserving reduced representation of the high-dimensional generalized coordinates and *(2)* the associated low-dimensional dynamic parameters. Specifically, our approach leverages a constrained AE (Otto et al., ) to learn a latent embedded submanifold in which physically-consistent low-dimensional dynamics are learned with a Lagrangian neural network (LNN) (Cranmer et al., ; Lutter & Peters, ).
Our approach differs from (Greydanus et al., ; Zhong & Leonard, ; Botev et al., ) in that it does not consider high-dimensional observations (images) of low-dimensional physical systems but systems with high-dimensional state spaces.
As *second* and *third contributions*, we reformulate both the LNN and the AE to account for the intrinsic geometry of their parameters. To this end, we introduce positive-definite layers in the LNN to parameterize the (reduced) mass-inertia matrix and use Riemannian optimization to infer the LNN parameters and biorthogonal AE weights.

Our approach infers physically-plausible models as both the ROM and the FOM conserve the reduced energy along reduced unforced Lagrangian trajectories and their embeddings. Moreover, it infers interpretable reduced-order dynamic parameters. We validate our approach by learning the dynamics of three simulated high-dimensional rigid and deformable systems: a pendulum, a rope, and a thin cloth.
Our results demonstrate that it efficiently learns reduced-order dynamics leading to accurate long-term predictions of high-dimensional systems.

## 2 Background

This section introduces the mathematical preliminaries and network architectures for learning reduced Lagrangian dynamics.
We adopt a differential-geometric perspective
to simultaneously investigate the use of nonlinear MOR projections and the structure preservation of Lagrangian systems.

### 2.1 Lagrangian Dynamics on the Configuration Manifold

We consider an $n$- degrees-of-freedom (DoF) mechanical system whose configuration space is identified with an $n$-dimensional smooth
manifold $\mathcal{Q}$ with a simple global chart.
Velocities $\dot{\bm{q}}$ at each configuration $\bm{q}\in\mathcal{Q}$ lie in the tangent space $\mathcal{T}_{\bm{q}}\mathcal{Q}$, an $n$-dimensional vector space composed of all vectors tangent to $\mathcal{Q}$ at $\bm{q}$.
The disjoint union of all tangent spaces $\mathcal{T}_{\bm{q}}\mathcal{Q}$ forms the tangent bundle $\mathcal{T}\mathcal{Q}$, a smooth $2n$-dimensional manifold.
The configuration manifold can be equipped with a Riemannian metric, i.e., a smoothly-varying
inner product acting on $\mathcal{T}\mathcal{Q}$. We consider the kinetic-energy metric, which, given a choice of local coordinates, equals the system’s mass-inertia matrix $\bm{M}(\bm{q})$.

A Lagrangian system is a tuple $(\mathcal{Q},\mathcal{L})$ of a Riemannian configuration manifold $\mathcal{Q}$ and a smooth time-independent Lagrangian function $\mathcal{L}:T\mathcal{Q}\to\mathbb{R}$.
The Lagrangian function is given by the difference between the system’s kinetic $T(\bm{q},\dot{\bm{q}})$ and potential $V(\bm{q})$ energies as

$$ $\mathcal{L}(\bm{q},\dot{\bm{q}})=T(\bm{q},\dot{\bm{q}})-V(\bm{q})=\frac{1}{2} \dot{\bm{q}}^{\intercal}\bm{M}(\bm{q})\dot{\bm{q}}-V(\bm{q}).$ (1) $$

Following the principle of least action, the equations of motion are given by the Euler-Lagrange equations
$\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\bm{q}}\right)-\frac{
\partial\mathcal{L}}{\partial\bm{q}}=\bm{\tau}$
with generalized forces $\bm{\tau}$.
These equations can be rewritten as

$$ $\bm{M}(\bm{q})\ddot{\bm{q}}+\bm{c}(\bm{q},\dot{\bm{q}})+\bm{g}(\bm{q})=\bm{ \tau},\quad\text{with}\quad\bm{c}(\bm{q},\dot{\bm{q}})=\left(\frac{d}{dt}\bm{M }(\bm{q})\right)\dot{\bm{q}}-\frac{1}{2}\left(\frac{\partial}{\partial\bm{q}} \left(\dot{\bm{q}}^{\intercal}\bm{M}(\bm{q})\dot{\bm{q}}\right)\right)^{ \intercal},$ (2) $$

and $\bm{g}(\bm{q})=\frac{\partial V(\bm{q})}{\partial\bm{q}}$, where $\bm{c}(\bm{q},\dot{\bm{q}})$ represents the influence of Coriolis forces.
Given a time interval $\mathcal{I}=\left[t_{0},t_{\text{f}}\right]$, trajectories $\bm{\gamma}:\mathcal{I}\to\mathcal{T}\mathcal{Q}:t\mapsto\left(\bm{q}(t),\dot{
\bm{q}}(t)\right)^{\mathsf{T}}$ of the Lagrangian system are obtained by solving the initial value problem (IVP)

$$ $\begin{cases}\frac{d}{dt}\bm{\gamma}\big{|}_{t}=\bm{X}\big{|}_{\bm{\gamma}(t)} =\left(\begin{smallmatrix}\dot{\bm{q}}(t)\\ \bm{M}^{-1}(\bm{q})\left(\bm{\tau}-\bm{c}(\bm{q},\dot{\bm{q}})-\bm{g}(\bm{q}) \right)\end{smallmatrix}\right)&\in\mathcal{T}_{\bm{\gamma}(t)}\mathcal{T} \mathcal{Q},,\[9.0pt] \bm{\gamma}(t_{0})=\left(\begin{smallmatrix}\bm{q}_{0}\\ \dot{\bm{q}}_{0}\end{smallmatrix}\right)&\in\mathcal{T}\mathcal{Q},\end{cases}$ (3) $$

with $\bm{X}\big{|}_{\bm{\gamma}(t)}$ the vector field defined by the Euler-Lagrange equations.
Note that, for Lagrangian systems, trajectories $\bm{\gamma}$ are lifted curves on the tangent bundle $\mathcal{T}\mathcal{Q}$.

### 2.2 Model Order Reduction

We consider a high-dimensional dynamical system, i.e., a FOM, which generates smooth trajectories $\bm{\gamma}:\mathcal{I}\to\mathcal{Q}$
characterized by an IVP of the form

$$ $\begin{cases}\frac{d}{dt}\bm{\gamma}\big{|}_{t}&=\bm{X}\big{|}_{\bm{\gamma}(t) }\in\mathcal{T}_{\bm{\gamma}(t)}\mathcal{Q},\quad t\in\mathcal{I},\[6.0pt] \bm{\gamma}(t_{0})&=\bm{\gamma}_{0}\in\mathcal{Q},\end{cases}$ (4) $$

where $\bm{X}\big{|}_{\bm{\gamma}(t)}$ is a smooth vector field defining the evolution of the system over time.
The goal of MOR is to accurately and efficiently approximate the set of solutions
$S=\left\{\bm{\gamma}(t)\in\mathcal{Q}\>|\>t\in\mathcal{I}\right\}\subseteq
\mathcal{Q}$
of the IVP (4).
To this end, MOR methods learn a reduced configuration manifold $\check{\mathcal{Q}}$ with $\text{dim}(\check{\mathcal{Q}})\!=\!d\ll\text{dim}(\mathcal{Q})\!=\!n$ and a ROM
$\frac{d}{dt}\check{\bm{\gamma}}\big{|}_{t}=\check{\bm{X}}\big{|}_{\check{\bm{
\gamma}}(t)}\in\mathcal{T}_{\check{\bm{\gamma}}(t)}\check{\mathcal{Q}}$.

Following Buchfink et al. (), we identify $\check{\mathcal{Q}}$ with an embedded submanifold $\varphi(\check{\mathcal{Q}})\subseteq\mathcal{Q}$ via a smooth embedding ${\varphi:\check{\mathcal{Q}}\to\mathcal{Q}}$. The reduced initial value $\check{\bm{\gamma}}_{0}\!=\!\rho(\bm{\gamma}_{0})$ and vector field ${\check{\bm{X}}\big{|}_{\check{\bm{\gamma}}(t)}=d\rho\big{|}_{\check{\bm{
\gamma}}(t)}\bm{X}\big{|}_{\bm{\gamma}(t)}}$ are defined via the point and tangent reduction maps $\rho:\mathcal{Q}\to\check{\mathcal{Q}}$ and $d\rho|_{\bm{q}}:\mathcal{T}_{\bm{q}}\mathcal{Q}\to\mathcal{T}_{\rho(\bm{q})}
\check{\mathcal{Q}}$ associated with $\varphi$.
The reduction maps must satisfy the projection properties

$$ $\rho\circ\varphi=\text{id}_{\check{\mathcal{Q}}}\quad\text{ and }\quad d\rho|_ {\varphi(\check{\bm{q}})}\circ d\varphi|_{\check{\bm{q}}}=\text{id}_{\mathcal{ T}_{\check{\bm{q}}}\check{\mathcal{Q}}},\quad\forall\check{\bm{q}}\in\check{ \mathcal{Q}}.$ (5) $$

Trajectories of the original system are then obtained via the approximation $\bm{\gamma}(t)\approx\varphi(\check{\bm{\gamma}}(t))$.
In this paper, we leverage structure-preserving MOR methods to identify a reduced-order Lagrangian $(\check{\mathcal{Q}},\check{\mathcal{L}})$, and obtain trajectories $\check{\bm{\gamma}}$ by solving the corresponding equations of motion.

### 2.3 Autoencoders for Model Order Reduction

A key aspect of MOR is the construction of the embedding $\varphi$ and corresponding point reduction $\rho$. AEs are well suited for representing these nonlinear mappings due to their expressivity and scalability (Lee & Carlberg, ; Otto et al., ).
An AE consists of an encoder network ${\rho:\mathbb{R}^{n}\to\mathbb{R}^{d}}$ that maps a high-dimensional vector $\bm{x}\in\mathbb{R}^{n}$ to a latent representation $\bm{z}\in\mathbb{R}^{d}$ with $d\ll n$, and a decoder network $\varphi:\mathbb{R}^{d}\to\mathbb{R}^{n}$ that reconstructs an approximation of the original data. Given a dataset $\{\bm{x}_{i}\}_{i=1}^{N}$, the AE parameters $\bm{\Psi}$ are trained by minimizing the reconstruction error

$$ $\ell_{\text{AE}}(\bm{\Psi})=\frac{1}{N}\sum_{i=1}^{N}\|\varphi\circ\rho(\bm{x} _{i})-\bm{x}_{i}\|^{2}.$ (6) $$

Notice that the projection properties (5) hold approximately when the loss (6) is small. However, in this paper, we leverage a constrained AE (Otto et al., ) that guarantees (5), which is essential when learning reduced representations of high-dimensional Lagrangian systems.

### 2.4 Lagrangian Neural Networks

The key idea of LNNs is to include physics-based inductive bias into deep networks, ensuring the learned dynamics models conserve energy and lead to physically-plausible trajectories. In this context, DeLaN (Lutter et al., ; Lutter & Peters, ) proposed to model the kinetic energy (via the mass-inertia matrix) and the potential energy of a Lagrangian function $\mathcal{L}$ as two networks $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})$ and $V(\bm{q};\bm{\theta}_{\text{V}})$ with parameters $\bm{\theta}=\{\bm{\theta}_{\text{T}},\bm{\theta}_{\text{V}}\}$. Trajectories of the learned dynamical system are then obtained via the equations of motion (2) by solving the IVP (3).
Given a set of $N$ observations $\left\{\bm{q}_{i},\dot{\bm{q}}_{i},\ddot{\bm{q}}_{i},\bm{\tau}_{i}\right\}_{i=
1}^{N}$, the networks are trained to minimize the loss

$$ $\ell_{\text{LNN}}(\bm{\theta})=\frac{1}{N}\sum_{i=1}^{N}\|\bm{f}(\bm{q}_{i}, \dot{\bm{q}}_{i},\bm{\tau}_{i};\bm{\theta})-\ddot{\bm{q}}_{i})\|^{2}+\lambda|| \bm{\theta}||_{2}^{2},$ (7) $$

where
${\bm{f}(\bm{q},\dot{\bm{q}},\bm{\tau};\bm{\theta})=\bm{M}^{-1}(\bm{q};\bm{
\theta}_{\text{T}})(\bm{\tau}-\bm{c}(\bm{q},\dot{\bm{q}})-\bm{g}(\bm{q};\bm{
\theta}_{\text{V}}))}$ and $\lambda$ is a regularization constant.
Following this approach, the total energy
$\mathcal{E}=T(\bm{q},\dot{\bm{q}})+V(\bm{q})=\frac{1}{2}\dot{\bm{q}}^{
\intercal}\bm{M}(\bm{q})\dot{\bm{q}}+V(\bm{q}),$
of unforced systems is conserved.
Next, we propose an enhanced geometric version of DeLaN, which we then use in Section 4 to learn the dynamics of a Lagrangian ROM.

## 3 Geometric Lagrangian Neural Networks

The mass-inertia matrix $\bm{M}(\bm{q})$ acts as the Riemannian metric on the configuration manifold $\mathcal{Q}$ for Lagrangian systems with quadratic kinetic energy structure. As such, $\bm{M}(\bm{q})$ belongs to the manifold of symmetric positive-definite (SPD) matrices $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$  (Bhatia, ; Pennec et al., ).
Existing LNNs enforce the symmetric positive-definiteness by predicting the Cholesky decomposition $\bm{L}(\bm{q})$ with a Euclidean network and describing the mass-inertia matrix as $\bm{M}=\bm{L}\bm{L}^{\mathsf{T}}$.
However, this parametrization overlooks the Riemannian geometry of both Cholesky and SPD spaces.
Here, we propose to learn the mass-inertia matrix via SPD networks, whose building blocks are formulated based on the Riemannian geometry of the SPD manifold.

### 3.1 Learning Positive-Definite Mass-Inertia Matrices on the SPD Manifold

Our goal is to learn a mapping $\bm{q}\mapsto\bm{M}(\bm{q}):\mathcal{Q}\to\mathcal{S}_{{\scriptscriptstyle{
\mathcal{++}}}}^{n}$, where $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}=\left\{\bm{\Sigma}\in
\text{Sym}^{n}\>\rvert\>\bm{\Sigma}\succ\mathbf{0}\right\}$ is the Riemannian manifold of SPD matrices and $\text{Sym}^{n}$ is the space of symmetric matrices (see App. A.1 for a background on the SPD manifold and associated Riemannian operations).
To this end, we build on recent works that generalize classical neural networks to operate directly on the SPD manifold by leveraging its gyrovectorspace structure (López et al., ; Nguyen, ; Nguyen et al., ). However, these networks consider that both inputs and outputs belong to $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$.

We propose a novel SPD network $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})=(g_{\mathcal{S}_{{\scriptscriptstyle{
\mathcal{++}}}}^{n}}\circ g_{\text{Exp}}\circ g_{\mathbb{R}})(\bm{q})$ with three consecutive components, as illustrated in Fig. 1-*left*.
First, the input configuration $\bm{q}\in\mathcal{Q}$ is fed to a standard Euclidean multilayer perceptron (MLP) $g_{\mathbb{R}}:\mathbb{R}^{n}\to\mathbb{R}^{n(n+1)/2}$, whose output elements are identified with those of a symmetric matrix $\bm{U}\in\text{Sym}^{n}$. Second, $\bm{U}$ is interpreted as an element of the tangent space $\mathcal{T}_{\bm{P}}\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$, and mapped onto the SPD manifold via an exponential map layer $\bm{X}=\text{Exp}_{\bm{P}}(\bm{U})$. Note that this type of layers is commonly used in hyperbolic neural networks (Shimizu et al., ). Third, $\bm{X}$ is fed to a set of $L_{\text{spd}}$ layers $g_{\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}^{(l)}:\mathcal{S}_{{
\scriptscriptstyle{\mathcal{++}}}}^{n}\to\mathcal{S}_{{\scriptscriptstyle{
\mathcal{++}}}}^{n}$, which outputs the mass-inertia matrix $\bm{M}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$.
We consider and compare two different fully-connected (FC) SPD layers — which are analogous to the FC layers in MLPs —, namely the gyrocalculus-based (GyroAI) layers of (López et al., ; Nguyen, ), and the gyrospace hyperplane-based (GyroSpd++) layers of (Nguyen et al., ). Both can be integrated with the ReEig layers (Huang & Gool, ) which are the SPD counterparts to ReLU-activations (see App. B for details on the SPD layers).
In total, the parameters $\bm{\theta}_{\text{T}}$ of our SPD network consist of the MLP parameters $\bm{\theta}_{\text{T},\mathbb{R}}$, and the SPD layer parameters $\bm{\theta}_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}\in
\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}\times\ldots\times\mathcal
{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$, including the basepoint $\bm{P}$ of the exponential map.

### 3.2 Model Training and Parameter Optimization

The proposed geometric LNN models the kinetic energy via the mass-inertia matrix parametrized as a SPD network $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})$ and the potential energy via a standard Euclidean MLP $V(\bm{q};\bm{\theta}_{\text{V}})$.
The network parameters $\bm{\theta}=\{\bm{\theta}_{\text{T}},\bm{\theta}_{\text{V}}\}$ are obtained as for classical LNNs by minimizing the loss (7).
However, as several parameters of the SPD network $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})$ belong to $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$, their geometry must be taken into account during training. To do so, we leverage Riemannian optimization methods (Absil et al., ; Boumal, ) to solve a problem of the form $\min_{\bm{x}\in\mathcal{M}}\ell(\bm{x})$ on a Riemannian manifold $\mathcal{M}$.
Conceptually, each iteration step in a first-order (stochastic) Riemannian optimization method consists of the three successive operations

$$ $\bm{\eta}_{t}\leftarrow h\big{(}\text{grad}\>\ell(\bm{x}_{t}),\bm{\tau}_{t-1} \big{)},\quad\bm{x}_{t+1}\leftarrow\text{Exp}_{\bm{x}_{t}}(-\alpha_{t}\bm{\eta }_{t}),\quad\bm{\tau}_{t}\leftarrow\mathrm{PT}_{\bm{x}_{t}\rightarrow\bm{x}_{t +1}}\big{(}\bm{\eta}_{t}\big{)}.$ (8) $$

First, the update $\bm{\eta}_{t}\in\mathcal{T}_{\bm{x}_{t}}\mathcal{M}$ is computed as a function $h$ of the Riemannian gradient $\text{grad}\>\ell$ at the estimate $\bm{x}_{t}$ and of $\bm{\tau}_{t-1}$, a parallel-transported previous update $\bm{\eta}_{t-1}\in\mathcal{T}_{\bm{x}_{t-1}}\mathcal{M}$ to the new $\mathcal{T}_{\bm{x}_{t}}\mathcal{M}$. Then, an update of the estimate $\bm{x}_{t}$ is obtained via projecting the update $\bm{\eta}_{t}$, that is scaled by a learning rate learning rate $\alpha_{t}$ onto the manifold with the exponential map.
Finally, the current update is parallel-transported to the tangent space of the updated estimate to prepare for the next iteration.
We use the Riemannian Adam (Becigneul & Ganea, ) implemented in Geoopt (Kochurov et al., ) to optimize the geometric LNN parameters. Notice that $\mathcal{M}$ is defined as a product of Euclidean and SPD manifolds to optimize the parameters $\bm{\theta}$.

Figure: Figure 1: Flowchart of the forward dynamics of the proposed reduced-order LNN. The reduction mappings and embeddings of the Lagrangian ROM are depicted in blue and parametrized via a constrained AE with biorthogonal layers (*right*). The ROM dynamics are learned via a latent geometric LNN depicted in orange. The mass-inertia matrix is parametrized via a SPD network (*left*).

## 4 Learning Reduced-order Lagrangian Dynamics

Like other gray-box models, our geometric LNN faces limitations in scaling to high-dimensional systems due to the increasing complexity of their dynamics.
In many cases, the solutions of high-dimensional equations of motion can be well approximated by a substantially lower-dimensional surrogate dynamic model.
Building on this assumption, we propose to reduce the problem dimensionality by learning a latent dynamical system that preserves the original Lagrangian structure.
The latent space and the reduced dynamics parameters are inferred jointly via structure-preserving MOR and a latent geometric LNN, as explained next. The complete forward dynamics of the proposed reduced-order LNN (RO-LNN) are illustrated in Fig. 1. Notice that we consider physical systems with high number of DoFs with observations $\{\bm{q}_{i},\dot{\bm{q}}_{i},\ddot{\bm{q}}_{i},\bm{\tau}_{i}\}_{i=1}^{T}$, and not high-dimensional representations of inherently low-dimensional state spaces as would be the case with images.

### 4.1 Lagrangian Reduced-order Model

Preserving the properties of the original FOM is crucial to ensure that the learned ROM displays similar behaviors. Therefore, we leverage the geometric framework introduced in Buchfink et al. () and learn a reduced Lagrangian $(\check{\mathcal{Q}},\check{\mathcal{L})}$ via structure-preserving MOR.
As Lagrangian trajectories are lifted curves $\bm{\gamma}(t)=\left(\bm{q}(t),\dot{\bm{q}}(t)\right)^{\mathsf{T}}$, the manifold to be reduced is the tangent bundle $\mathcal{T}\mathcal{Q}$. We define the lifted embedding ${\varphi:\mathcal{T}\check{\mathcal{Q}}\to\mathcal{T}\mathcal{Q}}$ for a smooth embedding $\varphi_{\mathcal{Q}}:\check{\mathcal{Q}}\to\mathcal{Q}$ as the pair

$$ $\varphi(\check{\bm{q}},\dot{\check{\bm{q}}})=(\varphi_{\mathcal{Q}}(\check{\bm {q}}),\;d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}\dot{\check{\bm{q}}}),$ (9) $$

with $d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}:\mathcal{T}_{\check{\bm{q}}}\check{
\mathcal{Q}}\to\mathcal{T}_{\varphi_{\mathcal{Q}}(\check{\bm{q}})}\mathcal{Q}$ denoting the pushforward, or differential, of $\varphi_{\mathcal{Q}}$ at $\check{\bm{q}}$. Analogously, we define a reduction map $\rho(\bm{q},\dot{\bm{q}})$ as the pair $(\rho_{\mathcal{Q}}(\bm{q}),d\rho_{\mathcal{Q}}|_{\bm{q}}\dot{\bm{q}})$ of point and tangent reductions $\rho_{\mathcal{Q}}:\mathcal{Q}\to\check{\mathcal{Q}}$ and $d\rho_{\mathcal{Q}}|_{\bm{q}}:\mathcal{T}_{\bm{q}}\mathcal{Q}\to\mathcal{T}_{
\rho_{\mathcal{Q}}(\bm{q})}\check{\mathcal{Q}}$.

Given $\varphi$ and $\rho$, the observed states $\{\bm{q}_{i},\dot{\bm{q}}_{i}\}_{i=1}^{T}$ are first mapped to the ROM via the point reduction mapping to obtain reduced initial values.
Then, the reduced Lagrangian function is constructed via the pullback of the lifted embedding as

$$ $\check{\mathcal{L}}=\varphi^{*}\mathcal{L}=\mathcal{L}\circ\varphi.$ (10) $$

The Euler-Lagrange equations of the reduced Lagrangian yield

$$ $\check{\bm{M}}(\check{\bm{q}})\ddot{\check{\bm{q}}}+\check{\bm{c}}(\check{\bm{ q}},\dot{\check{\bm{q}}})+\check{\bm{g}}(\check{\bm{q}})=\check{\bm{\tau}},$ (11) $$

where $\check{\bm{M}}(\check{\bm{q}})=d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}^{
\intercal}\;\bm{M}(\bm{q})\;d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}$, $\check{\bm{g}}(\check{\bm{q}})=d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}^{
\intercal}\;\bm{g}(\bm{q})$, $\check{\bm{\tau}}=d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}^{\intercal}\;\bm{\tau}$, and $\check{\bm{c}}(\check{\bm{q}},\dot{\check{\bm{q}}})$ is computed as function of $\check{\bm{M}}(\check{\bm{q}})$ as in (2).
The ROM (11) can be used to efficiently compute reduced trajectories $\check{\bm{\gamma}}(t)$ as solutions $\bm{\gamma}(t)\approx\varphi(\check{\bm{\gamma}}(t))$ of the FOM (2).

It is worth emphasizing that the preservation of the Lagrangian structure in the ROM guarantees the conservation of the reduced total energy $\check{\mathcal{E}}$ along the solutions $\check{\bm{\gamma}}(t)$ of (11) and the corresponding image curves $\varphi(\check{\bm{\gamma}}(t))$ as $\check{\mathcal{E}}=\mathcal{E}\circ\varphi$. We take advantage of this property and use the geometric LNN
to learn physically-consistent reduced dynamic parameters.
Specifically, the reduced mass-inertia matrix and the reduced potential energy of $\check{\mathcal{L}}$ are parametrized as two networks $\check{\bm{M}}(\check{\bm{q}};\bm{\theta}_{\check{\text{T}}})$ and $\check{V}(\check{\bm{q}};\bm{\theta}_{\check{\text{V}}})$ with parameters $\bm{\theta}=\{\bm{\theta}_{\check{\text{T}}},\bm{\theta}_{\check{\text{V}}}\}$.
Next, we discuss how to learn $\varphi_{\mathcal{Q}}$ and $\rho_{\mathcal{Q}}$.

### 4.2 Learning the Embedding and Point Reduction

For increased expressivity, we parametrize the point reduction $\rho_{\mathcal{Q}}$ and embedding $\varphi_{\mathcal{Q}}$ of our reduced Lagrangian as the encoder and decoder of an AE. To ensure that the reduction map $\rho(\bm{q},\dot{\bm{q}})$ satisfies the projection properties (5), we leverage the constrained AE architecture introduced by Otto et al. (). The encoder and decoder networks are given as a composition of feedforward layers ${\rho_{\mathcal{Q}}=\rho_{\mathcal{Q}}^{(1)}\circ\ldots\circ\rho_{\mathcal{Q}}
^{(L)}}$ and $\varphi_{\mathcal{Q}}=\varphi_{\mathcal{Q}}^{(L)}\circ\ldots\circ\varphi_{
\mathcal{Q}}^{(1)}$ with $\rho_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l}}\to\mathbb{R}^{n_{l-1}}$, $\varphi_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l-1}}\to\mathbb{R}^{n_{l}}$, and $d=n_{0}\leq\ldots\leq n_{L}=n$, where $d$ denotes the latent space dimension.
Notice that the pushforward $d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}$ and tangent reduction $d\rho_{\mathcal{Q}}|_{\bm{q}}$ are the differentials of the encoder and decoder networks.

The key to guarantee the projection properties (5) lies in the specific construction of layer pairs as

$$ $\rho_{\mathcal{Q}}^{(l)}(\bm{q}^{(l)})=\sigma_{-}\left(\bm{\Psi}_{l}^{\mathsf{ T}}(\bm{q}^{(l)}-\bm{b}_{l})\right)\quad\quad\text{and}\quad\quad\varphi_{ \mathcal{Q}}^{(l)}(\check{\bm{q}}^{(l-1)})=\bm{\Phi}_{l}\sigma_{+}(\check{\bm{ q}}^{(l-1)}+\bm{b}_{l}),$ (12) $$

where $(\bm{\Phi}_{l},\bm{\Psi}_{l})$ and $(\sigma_{+},\sigma_{-})$ are pairs of weight matrices and smooth activation functions, respectively, and $\bm{b}_{l}$ are bias vectors.
By constraining the pairs to satisfy $\bm{\Psi}_{l}^{\mathsf{T}}\bm{\Phi}_{l}=\mathbf{I}_{d}$ and $\sigma_{-}\circ\sigma_{+}=\text{id}$, each layer satisfies ${\rho_{\mathcal{Q}}^{(l)}\circ\varphi_{\mathcal{Q}}^{(l)}=\text{id}_{\mathbb{R
}^{n_{l-1}}}}$ and the constrained AE fullfills (5).
Otto et al. () adhere to the first constraint by defining pairs of biorthogonal matrices $(\bm{\Phi}_{l},\bm{\Psi}_{l})$ via an over-parametrization resulting in additional loss terms. Instead, we adopt a geometric approach that accounts for the Riemannian geometry of biorthogonal matrices. Specifically, we minimize the AE reconstruction error (6) via Riemannian optimization (see Equation (8)) by considering each pair $(\bm{\Phi}_{l},\bm{\Psi}_{l})$ as elements of the biorthogonal manifold $\mathcal{B}_{n_{l},n_{l-1}}=\{(\bm{\Phi},\bm{\Psi})\in\mathbb{R}^{n_{l}\times n
_{l-1}}\times\mathbb{R}^{n_{l}\times n_{l-1}}:\bm{\Psi}^{\mathsf{T}}\bm{\Phi}=
\bm{I}_{n_{l-1}}\}$ (see App. A.2 for a background on the biorthogonal manifolds and associated operations). As it will be shown later in our experiments, optimizing the biorthogonal weights on the biorthogonal manifold is crucial when jointly optimizing the latent space and the associated reduced-order dynamic parameters.
The second constraint is met by utilizing the smooth, invertible activation functions defined in (Otto et al., , Equation 12).
Additional details on the AE architecture, including activation functions, and layer derivatives, are provided in App. C.

### 4.3 Model Training

We finally propose to jointly learn the parameters $\bm{\Xi}\!=\!\{\bm{\Phi}_{l},\bm{\Psi}_{l},\bm{b}_{l}\}_{l=1}^{L}$ of the AE and $\bm{\theta}\!=\!\{\bm{\theta}_{\check{\text{T}},\mathbb{R}},\bm{\theta}_{
\check{\text{T}},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{d}},\bm{
\theta}_{\check{\text{V}}}\}$ of the latent geometric LNN. We consider two losses, both of which minimize the AE reconstruction error, the latent LNN loss, and a joint error on the reconstructed predictions.

Training on acceleration. Given a training set $\left\{\bm{q}_{i},\dot{\bm{q}}_{i},\ddot{\bm{q}}_{i},\bm{\tau}_{i}\right\}_{i=
1}^{N}$, the acceleration loss is

$$ $\ell_{\text{ROM, acc}}=\frac{1}{N}\sum_{i=1}^{N}\underbrace{\|\tilde{\bm{q}}_{ i}-\bm{q}_{i}\|^{2}+\|\tilde{\dot{\bm{q}}}_{i}-\dot{\bm{q}}_{i}\|^{2}+\|\tilde {\ddot{\bm{q}}}_{i}-\ddot{\bm{q}}_{i}\|^{2}}_{\ell_{\text{AE}}}+\underbrace{\| \ddot{\check{\bm{q}}}_{\text{p},i}-\ddot{\check{\bm{q}}}_{i}\|^{2}}_{\ell_{ \text{LNN},d}}+\underbrace{\|\tilde{\ddot{\bm{q}}}_{\text{p},i}-\ddot{\bm{q}}_ {i}\|^{2}}_{\ell_{\text{LNN},n}}+\lambda\|\bm{\theta}\|_{2}^{2}.$ (13) $$

$\tilde{\bm{q}}_{i}$, $\tilde{\dot{\bm{q}}}_{i}$, $\tilde{\ddot{\bm{q}}}_{i}$ denote the reconstructed position, velocity, and acceleration obtained by successively applying the lifted point reduction and embedding, and its derivatives.
The latent acceleration predictions $\ddot{\check{\bm{q}}}_{i}$ are predicted from the reduced equations of motion (11) as

$$ $\ddot{\check{\bm{q}}}_{\text{p}}={\bm{f}(\check{\bm{q}},\dot{\check{\bm{q}}}, \check{\bm{\tau}};\bm{\theta})=\check{\bm{M}}^{-1}(\check{\bm{q}};\bm{\theta}_ {\check{\text{T}}})(\check{\bm{\tau}}-\check{\bm{c}}(\check{\bm{q}},\dot{ \check{\bm{q}}})-\check{\bm{g}}(\check{\bm{q}};\bm{\theta}_{\check{\text{V}}}) )},$ (14) $$

The acceleration of the FOM is then computed as $\tilde{\ddot{\bm{q}}}_{\text{p}}=d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}\ddot
{\check{\bm{q}}}_{\text{p}}+d^{2}\varphi_{\mathcal{Q}}|_{(\check{\bm{q}},\dot{
\check{\bm{q}}})}\dot{\check{\bm{q}}}$, as shown in Fig. 7.

Training with multi-step integration.
The loss (13) requires the expensive computation of the second derivatives of the embedding. Moreover, it considers only single steps, while the learned dynamics are expected to predict multiple steps. Therefore, we also consider a multi-step loss that
which numerically integrates the latent acceleration predictions (14) before decoding. This is achieved achieved via $H$ Euler forward steps with constant integration time $\Delta t$.
Given sets of observations $\left\{\bm{q}_{i}(\mathcal{I}_{i}),\dot{\bm{q}}_{i}(\mathcal{I}_{i}),\bm{\tau}
_{i}(\mathcal{I}_{i})\right\}_{i=1}^{N}$ over intervals $\mathcal{I}_{i}=\left[t_{i},t_{i}+H\Delta t\right]$, the multi-step loss is

$$ $\begin{split}\ell_{\text{ROM,ODE}}=\frac{1}{HN}\sum_{i=1}^{N}\sum_{j=1}^{H}& \underbrace{\|\tilde{\bm{q}}_{i}(t_{i,j})-\bm{q}_{i}(t_{i,j})\|^{2}+\|\tilde{ \dot{\bm{q}}}_{i}(t_{i,j})-\dot{\bm{q}}_{i}(t_{i,j})\|^{2}}_{\ell_{\text{AE}}} \\ &+\underbrace{\|\dot{\check{\bm{q}}}_{\text{p},i}(t_{i,j})-\dot{\check{\bm{q}} }_{i}(t_{i,j})\|^{2}}_{\ell_{\text{LNN},d}}+\underbrace{\|\tilde{\dot{\bm{q}}} _{\text{p},i}(t_{i,j})-\dot{\bm{q}}_{i}(t_{i,j})\|^{2}}_{\ell_{\text{LNN},n}}+ \gamma\;\|\bm{\theta}\|_{2}^{2},\end{split}$ (15) $$

with latent velocity predictions $\dot{\check{\bm{q}}}_{\text{p},i}(t_{i,j})=\int_{t_{i}}^{t_{i,j}}\bm{f}\left(
\check{\bm{q}}_{i}(t),\dot{\check{\bm{q}}}_{i}(t),\check{\bm{\tau}}_{i}(t);\bm
{\theta}\right)$, velocity reconstructions $\tilde{\dot{\bm{q}}}_{\text{p},i}(t_{i,j})=d\varphi_{\mathcal{Q}}\left(\check{
\bm{q}}_{i}(t_{i,j})\right)$, and $t_{i,j}=t_{i}+j\Delta t$. Fig. 1 illustrates the resulting forward model.

## 5 Experiments

We first evaluate the geometric LNN on a simulated $2$-DoF planar pendulum. Second, we evaluate our geometric RO-LNN to learn the dynamics of three simulated high-dimensional systems: a $16$-DoF pendulum, a $192$-DoF rope, and a $600$-DoF thin cloth. Our experiments demonstrate the ability of our approach to learn reduced-order dynamics, resulting in accurate long-term predictions. Moreover, they highlight the importance of geometry as additional inductive bias in both LNN and AE. Details about simulation environments, datasets, network architectures, and training are provided in App. E for each experiment. Additional results are provided in App. F.

### 5.1 Learning Lagrangian Dynamics with Geometric LNNs

We start by evaluating the long-term prediction accuracy of the geometric LNN introduced in Section 3 on a $2$-DoF planar pendulum. Train and test trajectories are obtained by solving its equations of motion with randomly sampled initial positions and zero velocities.
We compare the performance of our geometric LNN against DeLaN, which parametrizes the mass-inertia matrix via a Cholesky network. We consider two variants, where the kinetic and potential energy networks are independent (Lutter & Peters, ) or share parameters ($\bm{\theta}_{\text{T}}\cap\bm{\theta}_{\text{V}}$) (Lutter et al., ). Moreover, we evaluate the performance of different architectures for our SPD network $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})$. First, we compare two exponential map layers, where the basepoint $\bm{P}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ is either set as $\mathbf{I}$ or learned as a parameter. Second, we evaluate the influence of the SPD layers by comparing performances without any SPD layers, and with GyroAI+ReEig, GyroSpd++, or GyroSpd+++ReEig layers.

Fig. 2-*left* shows the acceleration prediction errors on testing data for selected architectures trained on acceleration data.
We observe that our geometric LNNs consistently outperform both DeLaNs, especially in low-data regimes.
Interestingly, the geometric LNNs with SPD layers do not noticeably outperform those employing solely Euclidean and exponential-map layers (see App. F.1.3 for an extensive analysis on SPD layers and App. F.1.5 for training times).
Fig. 2-*middle, right* depict long-term trajectory predictions obtained via Euler integration of the predicted state derivatives. In addition to the aforementioned models, we consider several geometric LNNs trained with multi-step integration.
We observe that the geometric LNNs lead to significantly better long-term predictions than the DeLaNs, with the SPD layers following the previously-observed performance trend (see also App. F.1.1). Moreover, the geometric LNN trained over $8$ integration steps outperforms those trained on accelerations, and results in the most accurate predictions (see also App. F.1.4).
Overall, our results validate the effectiveness of considering additional inductive bias given by *(1)* the intrinsic geometry of the mass-inertia matrix, and *(2)* multi-step integration when training LNNs. These results are further validated on a second dataset in App. F.1.2.

Figure: Figure 2: $2$-DoF pendulum: *Left*: Median acceleration prediction error for different LNNs and training set sizes $\mathcal{D}_{\text{train}}$ over $10$ test trajectories. Shaded regions represent first and third quartiles. *Middle, right:* Numerically-integrated position predictions and ground truth trajectory ( ).

### 5.2 Learning Reduced-Order Lagrangian Dynamics

Next, we learn reduced-order dynamics of several simulated high-dimensional physical systems.
Due to the limited performance improvements and increased computational complexity of the SPD layers, we only evaluate SPD networks $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})=(g_{\text{Exp}}\circ g_{\mathbb{R}})(\bm
{q})$ with $L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}=0$.

#### 5.2.1 Coupled Pendulum ( 16 16 16 16 DoFs )

As a first ROM example, we consider a $16$-DoF pendulum.
We ensure that the dynamics are reducible by constraining the motion of the last $12$-DoF as nonlinear combinations of the first $4$-DoF. As such, the high-dimensional pendulum is reducible to $4$ latent dimensions.

Learning High-Dimensional Dynamics.
We model the high-dimensional dynamics with our RO-LNN. The parameters of the structure-preserving AE and latent geometric LNN are jointly trained via Riemannian optimization. We set the latent space dimension to $d=4$ and the exponential map basepoint to $\mathbf{I}$ (see App. E.2 for implementation details).
We consider models trained on acceleration data via (13), and with multi-step integration using (15) with $H_{\text{train}}=8$.
Fig. 3 shows the AE reconstruction and RO-LNN predictions ($H_{\text{test}}=25$) for selected DoFs of a test trajectory. Average testing errors for different $H_{\text{test}}$ are reported in App. F.2.1. We observe that the RO-LNN trained with multi-step integration consistently outperforms the RO-LNN trained on accelerations, validating the importance of considering successive states during training.
Moreover, we compare our RO-LNN with full-order LNNs that directly learn high-dimensional dynamic parameters. As shown in Table 1, even the best-performing FOM produces prediction errors that are orders of magnitude higher than the RO-LNN errors. Such erroneous acceleration predictions did not allow us to obtain stable velocity and position predictions.

Figure: Figure 3: $16$-DoF pendulum: Comparisons of the position and velocity predictions from the RO-LNNs trained on acceleration ( ) and via multi-step integration ( ) with the ground truth ( ). The corresponding AE reconstructions ( ) and ( ) are depicted for completeness. The model is updated with a new initial condition $(\bm{q}_{0},\dot{\bm{q}}_{0})^{\mathsf{T}}$ every $0.025\text{\,}\mathrm{s}$ ($H_{\text{test}}=25$).

**Table 1: Mean and standard deviation of the acceleration prediction errors over $10$ test trajectories.**
| LNN | RO-LNN (acc) | RO-LNN (ODE) |
| --- | --- | --- |
| $3.80e^{8}\pm 1.12e^{9}$ | $2.10e^{3}\pm 2.66e^{3}$ | $\bm{6.80e^{2}\pm 1.63e^{3}}$ |

Joint training.
It is important to emphasize that the quality of the AE-reconstructed states is crucial for learning accurate dynamics. If the ROM cannot effectively capture the solution space of the FOM, the learned dynamics may systematically
deviate from the ground truth.
Here, we analyze the influence of the AE on the overall performance of our RO-LNN. Specifically, we compare our model with the joint training process with *(1)* a ROM which sequentially trains the independent constrained AE and latent LNN, and *(2)* a jointly-trained ROM that employs the overparametrization of (Otto et al., ), and thus overlooks the geometry of the biorthogonal AE weights.
All models are trained with $d=4$ via multi-step integration with $H_{\text{train}}=8$ and $\Delta t=10^{-3}$\text{\,}\mathrm{s}$$.

Figure 4 shows AE reconstruction errors (position and velocity) along with the high-dimensional and latent prediction errors averaged over $10$ testing trajectories. We observe that the jointly-trained models outperform the sequentially-trained one. Although the latter learns accurate dynamics in the latent space, the reconstructed high-dimensional dynamics result in significantly higher prediction errors (see the right $y$-scale).
Moreover, the overparametrized model leads to higher errors in all loss components compared to the RO-LNN trained on the biorthogonal manifold. The overparametrized AE is also noticeably harder to train jointly with the latent LNN: In our experiments, $4$ out of $6$ trainings led to nan-values. Altogether, our results suggest that considering the intrinsic geometry of the biorthogonal manifold enhances the performance and stability of the parameter optimization (see App. F.2.2 for further comparisons).

Figure: Figure 4: Median and quartiles of the errors of jointly-trained (on-biorthogonal-manifold ( ) vs overparametrized ( )) and sequentially-trained ( ) models. Note that the latter have their own $y$-axis in the left three plots, which is at least an order of magnitude higher.

#### 5.2.2 Rope ( 192 192 192 192 DoFs )

Next, we consider high-dimensional deformable systems and learn the dynamics of a $192$-dimensional rope. Trajectories are generated by tracing a circle with one end while keeping the other fixed. We train our RO-LNN with $d=10$ via multi-step integration (see App. E.3 for details on the data and implementation).
Fig 5 shows the predicted rope configurations for a horizon of $0.025\text{\,}\mathrm{s}$ ($H_{\text{test}}=25$). Our model accurately predicts the high-dimensional dynamics of the rope (see also the predictions for selected DoFs in App. F.3).
Notice that our attempts to train a full-order LNN were hindered by the high dimensionality of the system.

Figure: Figure 5: Predicted rope position ( ) and ground truth ( ) at selected timesteps. The grey circle depicts the circular trajectory of the end of the rope.
Refer to caption: extracted/5951634/figures/rope/rope_frames_time_png.png

Latent space dimension.
Selecting an appropriate latent dimension is an important consideration for MOR.
Here, we study the influence of the latent dimension on the performance of our RO-LNN.
Table 2 reports average reconstruction and prediction errors of RO-LNNs with latent dimensions $d=\{4,6,10,14\}$ on $10$ testing trajectories.
We observe that errors initially decrease as the latent dimension increases, suggesting that higher-dimensional latent spaces better capture the original high-dimensional dynamics. However, the errors increase beyond a certain latent dimension, indicating that the latent LNN becomes harder to train. In other words, the choice of latent dimension trades off between latent space expressivity and the limitations of LNNs in higher dimensions.

**Table 2: Rope reconstruction errors and prediction loss components from (15) for RO-LNN with different latent dimensions.**
| $d$ | $\|\tilde{\bm{q}}-\bm{q}\|^{2}$ | $\|\tilde{\dot{\bm{q}}}-\dot{\bm{q}}\|^{2}$ | $\ell_{\text{LNN},n}$ | $\ell_{\text{LNN},d}$ |
| --- | --- | --- | --- | --- |
| $4$ | $2.95e^{-2}\pm 3.24e^{-2}$ | $4.63e^{-2}\pm 7.69e^{-2}$ | $5.06^{-2}\pm 8.23e^{-2}$ | $3.00e^{-3}\pm 1.34e^{-2}$ |
| $6$ | $8.19e^{-3}\pm 9.11e^{-3}$ | $1.52e^{-2}\pm 2.92e^{-2}$ | $2.40e^{-2}\pm 4.57e^{-2}$ | $\bm{1.81e^{-3}\pm 5.69e^{-3}}$ |
| $10$ | $\bm{4.75e^{-3}\pm 5.73e^{-3}}$ | $\bm{6.53e^{-3}\pm 1.62e^{-2}}$ | $\bm{1.43e^{-2}\pm 3.45e^{-2}}$ | $2.87e^{-3}\pm 8.92e^{-3}$ |
| $14$ | $6.16e^{-2}\pm 6.52e^{-2}$ | $5.66e^{-2}\pm 8.48e^{-2}$ | $9.22e^{-1}\pm 1.19e^{-1}$ | $5.06e^{-2}\pm 1.08e^{-1}$ |

#### 5.2.3 Cloth ( 600 600 600 600 DoFs )

As a final example, we learn the dynamics of an even higher-dimensional deformable system: A simulated $600$-DoF thin cloth falling onto a sphere. We generate different trajectories by randomly varying the sphere radius, which lead to different actuations applied to the deformable structure. We train our RO-LNN with $d=10$ via multi-step integration on a set of $|\mathcal{D_{\text{train}}}|=8000$ samples (see App. E.4 for details on data and implementation).
Fig. 6-*middle* depicts the predicted cloth configurations for a horizon of $0.025\text{\,}\mathrm{s}$ ($H_{\text{test}}=25$), showing that our model accurately predicts the cloth high-dimensional dynamics (see also the predictions for selected DoFs and errors in App. F.4).

Figure: Figure 6: Predicted cloth configuration at selected times for $25$-steps and $2500$-steps horizons.

Full-horizon predictions.
We evaluate prediction quality of our model over a longer time horizon. Namely, given an initial state $(\bm{q}_{0},\dot{\bm{q}}_{0})^{\mathsf{T}}$, we use our model to predict the motion of the cloth for $H_{\text{test}}=2500$ timesteps.
Fig. 6-*bottom* shows the cloth configurations predicted over the full horizon (see App. F.4 for predictions of selected DoFs and energy).
Our model accurately predicts the cloth high-dimensional dynamics over the full horizon, suggesting that it successfully learned the constraints and actuation effects of the cloth falling on a sphere.

## 6 Discussion

This paper proposed a novel geometric network architecture for learning the dynamics of high-dimensional Lagrangian systems.
By leveraging geometry and physics as inductive bias, our model is physically-consistent and conserves the reduced energy, infers interpretable reduced dynamic parameters, and effectively learns reduced-order dynamics that accurately describe the behavior of high-dimensional rigid and deformable systems.
It is worth emphasizing the ubiquitous role of geometry in our model. First, considering the geometry of Lagrangian systems allows us to define nonlinear embeddings that ensure the preservation of the energy along solution trajectories and the interpretability of the model. Second, considering the geometry of full-order LNNs, latent LNNs, and AE parameters during their (joint) optimization leads to increased performance, especially at low-data regimes, in addition to guaranteeing that they belong to the manifold of interest.

It is important to note that the prediction quality can only be as good as the AE reconstruction. In other words,
if the ROM cannot effectively capture the FOM solution space, the learned dynamics will systematically deviate from the ground truth.
In this sense, the selection of the latent dimension is crucial, but non-trivial due to the unknown degree of coupling in the FOM.
In future work, we plan to extend our approach to applications beyond system identification. In particular, we will work on the design control strategy based on the learned reduced dynamics for robotic manipulation of deformable objects and dynamic control of soft robots.

#### Acknowledgments

This work was supported by ERC AdV grant BIRD, Knut and Alice Wallenberg Foundation and Swedish Research Council.

## Appendix A Riemannian Manifolds

### A.1 The Manifold of SPD Matrices

In this section, we provide a brief overview on the Riemannian manifold of symmetric positive-definite (SPD) matrices.
The set of SPD matrices is defined as $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}=\left\{\bm{\Sigma}\in
\text{Sym}^{n}\>\rvert\>\bm{\Sigma}\succ\mathbf{0}\right\}$ with $\text{Sym}^{n}$ denotes the set of symmetric $\mathbb{R}^{n\times n}$ matrices. The set $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ forms a smooth manifold of dimension $\text{dim}(\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n})=\frac{n(n+1)}
{2}$, which can be represented as the interior of a convex cone embedded in $\text{Sym}^{n}$.
The tangent space $\mathcal{T}_{\bm{\Sigma}}\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ at any point ${\bm{\Sigma}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}$ is identified with the symmetric space $\text{Sym}^{n}$ with origin at $\bm{\Sigma}$.

The SPD manifold can be endowed with various Riemannian metrics. In this paper, we consider the affine-invariant metric (Pennec et al., ), which is widely used due to its excellent theoretical properties. The affine-invariant metric defines the inner product $g:\mathcal{T}_{\bm{\Sigma}}\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n
}\times\mathcal{T}_{\bm{\Sigma}}\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}
}}^{n}\to\mathbb{R}$ given, for two matrices $\bm{T}_{1}$, $\bm{T}_{2}\in\mathcal{T}_{\bm{\Sigma}}\mathcal{S}_{{\scriptscriptstyle{
\mathcal{++}}}}^{n}$, as

$$ $\langle\bm{T}_{1},\bm{T}_{2}\rangle_{\bm{\Sigma}}\;\;=\;\;\text{tr}(\bm{\Sigma }^{-\frac{1}{2}}\bm{T}_{1}\bm{\Sigma}^{-1}\bm{T}_{2}\bm{\Sigma}^{-\frac{1}{2}}).$ (16) $$

The affine-invariant metric defines the following geodesic distance

$$ $d_{\mathcal{M}}(\bm{\Lambda},\bm{\Sigma})=\|\log(\bm{\Sigma}^{-\frac{1}{2}}\bm {\Lambda}\bm{\Sigma}^{-\frac{1}{2}})\|_{\text{F}}.$ (17) $$

To map back and forth between the Euclidean tangent space and the manifold, we can employ the corresponding exponential and logarithmic maps $\text{Exp}_{\bm{\Sigma}}:\mathcal{T}_{\bm{\Sigma}}\mathcal{S}_{{
\scriptscriptstyle{\mathcal{++}}}}^{n}\to\mathcal{S}_{{\scriptscriptstyle{
\mathcal{++}}}}^{n}$ and $\text{Log}_{\bm{\Sigma}}:\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}
\to\mathcal{T}_{\bm{\Sigma}}\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$.
The exponential and logarithmic maps are then computed in closed form as

$$ $\displaystyle\bm{\Lambda}$ $\displaystyle=\text{Exp}_{\bm{\Sigma}}(\bm{L})=\bm{\Sigma}^{\frac{1}{2}}\exp( \bm{\Sigma}^{-\frac{1}{2}}\bm{L}\bm{\Sigma}^{-\frac{1}{2}})\bm{\Sigma}^{\frac{ 1}{2}},$ (18) $\displaystyle\bm{L}$ $\displaystyle=\text{Log}_{\bm{\Sigma}}(\bm{\Lambda})=\bm{\Sigma}^{\frac{1}{2}} \log(\bm{\Sigma}^{-\frac{1}{2}}\bm{\Lambda}\bm{\Sigma}^{-\frac{1}{2}})\bm{ \Sigma}^{\frac{1}{2}},$ (19) $$

where $\exp(\cdot)$ and $\log(\cdot)$ denote the matrix exponential and logarithm functions.
The parallel transport $\mathrm{PT}_{\bm{\Sigma}\rightarrow\bm{\Lambda}}:\mathcal{T}_{\bm{\Sigma}}
\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}\to\mathcal{T}_{\bm{
\Lambda}}\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ allowing us to move elements between different tangent spaces is defined as

$$ $\bm{\tilde{T}}=\mathrm{PT}_{\bm{\Sigma}\rightarrow\bm{\Lambda}}\big{(}\bm{T} \big{)}=\bm{A}_{\bm{\Sigma}\to\bm{\Lambda}}\;\bm{T}\;\bm{A}_{\bm{\Sigma}\to\bm {\Lambda}}^{\mathsf{T}},$ (20) $$

with $\bm{A}_{\bm{\Sigma}\to\bm{\Lambda}}=\bm{\Lambda}^{\frac{1}{2}}\bm{\Sigma}^{-
\frac{1}{2}}$.
The exponential map, logarithmic map, and parallel transport are key operations for optimizing parameters on the SPD manifold in Sections 3 and 4.

### A.2 The Biorthogonal Manifold

Pairs of matrices matrices $\bm{\Phi},\bm{\Psi}\in\mathbb{R}^{n\times d}$ with $n\geq d\geq 1$ that fullfil $\bm{\Psi}^{\mathsf{T}}\bm{\Phi}=\mathbf{I}$ form the smooth manifold $\mathcal{B}_{n,d}=\{(\bm{\Phi},\bm{\Psi})\in\mathbb{R}^{n\times d}\times
\mathbb{R}^{n\times d}\>\rvert\>\bm{\Psi}^{\mathsf{T}}\bm{\Phi}=\bm{I}_{d}\}$ (Otto et al., ). The biorthogonal manifold $\mathcal{B}_{n,d}$ is an embedded submanifold of the Euclidean product space $\mathbb{R}^{n\times d}\times\mathbb{R}^{n\times d}$ of dimension $\text{dim}\left(\mathcal{B}_{n,d}\right)=2nd-d^{2}$.
The tangent space at a point $\left(\bm{\Phi},\bm{\Psi}\right)\in\mathcal{B}_{n,d}$ is given by

$$ $\mathcal{T}_{\left(\bm{\Phi},\bm{\Psi}\right)}\mathcal{B}_{n,d}=\left\{\left( \bm{V},\bm{W}\right)\in\mathbb{R}^{n\times d}\times\mathbb{R}^{n\times d}\;:\; \bm{W}^{\mathsf{T}}\bm{\Phi}+\bm{\Psi}^{\mathsf{T}}\bm{V}=\mathbf{0}\right\}.$ (21) $$

Vector can be projected from the embedding space onto the tangent space via the projection operation $\text{Proj}_{\left(\bm{\Phi},\bm{\Psi}\right)}:\mathbb{R}^{n\times d}\times
\mathbb{R}^{n\times d}\to\mathcal{T}_{\left(\bm{\Phi},\bm{\Psi}\right)}
\mathcal{B}_{n,d}$ is given by

$$ $\text{Proj}_{\left(\bm{\Phi},\bm{\Psi}\right)}\left(\bm{X},\bm{Y}\right)=\left (\bm{X}-\bm{\Psi}\bm{A},\bm{Y}-\bm{\Phi}\bm{A}^{\mathsf{T}}\right),$ (22) $$

where $\bm{A}$ is the solution of the Sylvester equation
$\bm{A}(\bm{\Phi}^{\mathsf{T}}\bm{\Phi})+(\bm{\Psi}^{\mathsf{T}}\bm{\Psi})\bm{A
}=\bm{Y}^{\mathsf{T}}\bm{\Phi}+\bm{\Psi}^{\mathsf{T}}\bm{X}$.

Riemannian operations on the biorthogonal manifold are generally difficult to obtain. Therefore, we leverage a retraction $\mathrm{R}_{\left(\bm{\Phi},\bm{\Psi}\right)}:\mathcal{T}_{\left(\bm{\Phi},\bm
{\Psi}\right)}\mathcal{B}_{n,d}\to\mathcal{B}_{n,d}$ as first-order approximation for the exponential map. The retraction is formulated as

$$ $\mathrm{R}_{\left(\bm{\Phi},\bm{\Psi}\right)}\left(\bm{V},\bm{W}\right)=\left( (\bm{\Phi}+\bm{V})\left[(\bm{\Psi+\bm{W}})^{\mathsf{T}}(\bm{\Phi}+\bm{V}) \right]^{-1},\bm{\Psi}+\bm{W}\right).$ (23) $$

A first-order approximation of the parallel transport operation on the biorthogonal manifold $\mathrm{PT}_{\left(\bm{\Phi}_{1},\bm{\Psi}_{1}\right)\rightarrow\left(\bm{\Phi
}_{2},\bm{\Psi}_{2}\right)}:\mathcal{T}_{\left(\bm{\Phi}_{1},\bm{\Psi}_{1}
\right)}\mathcal{B}_{n,d}\to\mathcal{T}_{\left(\bm{\Phi}_{2},\bm{\Psi}_{2}
\right)}\mathcal{B}_{n,d}$ is then defined via the successive application of the retraction and projection as
$\mathrm{PT}_{\left(\bm{\Phi}_{1},\bm{\Psi}_{1}\right)\rightarrow\left(\bm{\Phi
}_{2},\bm{\Psi}_{2}\right)}=\text{Proj}_{\left(\bm{\Phi}_{2},\bm{\Psi}_{2}
\right)}\circ\mathrm{R}_{\left(\bm{\Phi}_{1},\bm{\Psi}_{1}\right)}$.

## Appendix B Layer types in SPD-Networks

Here, we briefly introduce the different layers used in the SPD networks of Section 3.

Euclidean Layers. Our SPD network leverages classical fully-connected layers to model functions that return elements on the tangent space of a manifold. The output of the $l$-th Euclidean layer $\bm{x}^{(l)}$ is given by

$$ $\bm{x}^{(l)}=\sigma(\bm{A}_{l}\bm{x}^{(l-1)}+\bm{b}_{l}),$ (24) $$

with $\bm{A}_{l}\in\mathbb{R}^{n_{l}\times n_{l-1}}$ and $\bm{b}_{l}\in\mathbb{R}^{n_{(l)}}$ the weight matrix and bias of the layer $l$, and $\sigma$ a nonlinear activation function of choice.

Exponential Map Layers. In the case of SPD manifolds, the exponential map layer is used to map layer inputs $\bm{X}^{(l-1)}\in\text{Sym}^{n}$ from the tangent space to the manifold. The layer output is given by

$$ $\bm{X}^{(l)}=\text{Exp}_{\bm{P}}(\bm{X}^{(l-1)}),$ (25) $$

with $\bm{P}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ denoting the basepoint of the considered tangent space. One option is to define $\bm{P}$ as equal to the identity matrix $\mathbf{I}$. In that case, the layer input is assumed to lie in the tangent space at the origin of the cone. As the tangent spaces only give a local approximation to the manifold’s curvature, the basepoint can instead be learned as a network parameter $\bm{P}\in\bm{\theta}_{\text{T, $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}
}}^{n}$}}$ to potentially increase expressivity.

Gyrocalculus-based FC SPD Layers (GyroAI). Together with gyrocalculus formulations for the SPD manifold, López et al. () and Nguyen () presented a version of fully connected layers for SPD data. With $\oplus$ as addition and $\otimes$ as scaling operations for SPD gyrovectors,

$$ $\bm{X}^{(l)}=\left(\bm{A}_{l}\otimes\bm{X}^{(l-1)}\right)\oplus\bm{B}_{l},$ (26) $$

defines the $l$-th gyrolayer output $\bm{X}^{(l)}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ for the input $\bm{X}^{(l-1)}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$. Here,
the gyrovector operations correspond to $\bm{A}\otimes\bm{X}=\text{exp}(\bm{A}*\text{log}(\bm{X}))$ with $*$ denoting pointwise multiplication and $\bm{X}\oplus\bm{B}=\bm{X}^{\frac{1}{2}}\bm{B}\bm{X}^{\frac{1}{2}}$. Matrices $\bm{A}_{l},\bm{B}_{l}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ are part of the network parameters $\bm{\theta_{\text{T, $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$}}}$.

Gyrospace hyperplane-based FC SPD Layers (GyroSpd++). Inspired by their equivalent in hyperbolic spaces (Shimizu et al., ), another version of a gyrocalculus-based FC layer for SPD networks was proposed by Nguyen et al. (). Here, the classic combination of matrix product with a weight and addition of a bias is treated as hyperplane equation. Therefore, these fully-connected layers compute
the distance to a hyperplane that is located in the origin of the manifold from a specific point $\bm{V}$, that is obtained from the layer input $\bm{X}^{(l-1)}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ via the gyrocalculus formulation of the hyperplane equation

$$ $\bm{V}(\bm{X}^{(l-1)})=\langle\ominus\bm{B}_{l}\oplus\bm{X}^{(l-1)},\bm{A}_{l}\rangle,$ (27) $$

that is parametrized by the network parameters $\bm{A}_{l},\bm{B}_{l}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$. The additive inverse $\ominus$ fulfills $\ominus\bm{X}=\bm{X}^{-1}$.
The final layer output is given by

$$ $\bm{X}^{(l)}=\text{exp}((\bm{X}^{(l-1)})*\bm{I}_{\text{v}}),$ (28) $$

where $\bm{I}_{\text{v}(i,j)}=\frac{1}{\sqrt{2}}$ for $i\neq j$ and $\bm{I}_{\text{v}(i,j)}=1$ for $i=j$.

ReEig Nonlinearities. The ReEig layers were first introduced by Huang & Gool (), and were extended to a more general context in Nguyen (). Even though $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ already classifies as nonlinear space, this type of layer adds the option of adding further nonlinearities to the network. The layer output is defined as

$$ $\bm{X}^{(l)}=\sigma_{\text{spd}}(\bm{X}^{(l-1)})=\bm{U}\text{diag}(\text{max}( \epsilon\mathbf{I},\bm{V}))\bm{U}^{\mathsf{T}},$ (29) $$

with eigenvalue decomposition of $\bm{X}^{(l-1)}=\bm{U}\text{diag}(\bm{V})\bm{U}^{\mathsf{T}}$. The constant $\epsilon$ acts as rectification threshold to avoid smaller eigenvalues and hereby regularize the output.
Under the affine-invariant distance, this nonlinearity layer is analogous to a ReLu nonlinearity in the Euclidean network layers.

## Appendix C Constrained AE: Implementation Details

The constrained AE is implemented following (Otto et al., ).
To guarantee the projection properties, the nonlinear activation functions employed in the encoder and decoder network $\sigma_{-}$ and $\sigma_{+}$ must satisfy $\sigma_{-}\circ\sigma_{+}=\text{id}$. To do so, they are defined by

$$ $\sigma_{\pm}(x_{i})=\frac{bx_{i}}{a}\mp\frac{\sqrt{2}}{a\sin(\alpha)}\pm\frac{ 1}{a}\sqrt{\left(\frac{2x_{i}}{\sin(\alpha)\cos(\alpha)}\mp\frac{\sqrt{2}}{ \cos(\alpha)}\right)+2a},$ (30) $$

with

$$ $\begin{cases}a&=\csc^{2}(\alpha)-\sec^{2}(\alpha),\\ b&=\csc^{2}(\alpha)+\sec^{2}(\alpha).\end{cases}$ (31) $$

The activations then resemble smooth, rotation-symmetric versions of the common leaky ReLu activations. $0<\alpha<\frac{\pi}{4}$ is a parameter to set the slope of the activation functions. Throughout our experiments, we set $\alpha=\frac{\pi}{8}$.

For each layer $l$, the weight matrices are given by the pair $\left(\bm{\Phi}_{l},\bm{\Psi}_{l}\right)\in\mathcal{B}_{n_{l},n_{l-1}}$ (see App. A.2). In the proposed RO-LNN, we train these weights via Riemannain optimization on the biorthogonal manifold, analogous to the Riemannian optimization on the SPD manifold presented in 3.2. Notice that the optimization of the AE parameters $\bm{\Xi}\!=\!\{\bm{\Phi}_{l},\bm{\Psi}_{l},\bm{b}_{l}\}_{l=1}^{L}$ takes place on a product of Euclidean and biorthogonal manifold.

When learning dynamic latent space for Lagrangian MOR, the AE should learn the lifted embedding (9) and corresponding reduction maps. Therefore, the outputs of the encoder and decoder networks are differentiated with respect to their inputs. In our implementation, we take layerwise analytical derivatives and obtain the full differentials via the chain rule. A second derivative of the encoder and decoder networks enables training on acceleration data. Full-dimensional acceleration data is encoded via the reduction map $\rho_{\mathcal{Q}}$ as $\ddot{\check{\bm{q}}}=d\rho_{\mathcal{Q}}|_{\bm{q}}\ddot{\bm{q}}+d^{2}\rho_{
\mathcal{Q}}|_{(\bm{q},\dot{\bm{q}})}\dot{\bm{q}}$, and latent accelerations are reconstructed via the embedding $\varphi_{\mathcal{Q}}$ as $\tilde{\ddot{\bm{q}}}=d\varphi_{\mathcal{Q}}|_{\check{\bm{q}}}\ddot{\check{\bm
{q}}}+d^{2}\varphi_{\mathcal{Q}}|_{(\check{\bm{q}},\dot{\check{\bm{q}}})}\dot{
\check{\bm{q}}}$.
This requires the expensive computation of second-order derivatives, but is nevertheless part of our experiments.

Given a dataset of $T$ observations $\left\{\bm{q}_{i},\dot{\bm{q}}_{i},\ddot{\bm{q}}_{i}\right\}_{i=1}^{T}$, the AE is trained on a loss

$$ $\ell_{\text{cAE}}=\frac{1}{N}\sum_{i=1}^{N}\|\tilde{\bm{q}}_{i}-\bm{q}_{i}\|^{ 2}+\|\tilde{\dot{\bm{q}}}_{i}-\dot{\bm{q}}_{i}\|^{2}+\|\tilde{\ddot{\bm{q}}}_{ i}-\ddot{\bm{q}}_{i}\|^{2},$ (32) $$

of position reconstructions $\tilde{\bm{q}}$, velocity reconstructions $\tilde{\dot{\bm{q}}}$, and acceleration reconstruction $\tilde{\ddot{\bm{q}}}$. If the velocity or acceleration encodings are not relevant, the respective terms can simply be removed.

We compare the performance of the constrained AE trained via Riemannian optimization with the training scheme introduced in (Otto et al., ), where the biorthogonality of weights is achieved via the introduction of additional penalty losses. In this case, the overall loss only reaches a minimum value if the overparametrized weight matrices actually fulfill the biorthogonality conditions for each layer $l$.
The overparametrized constrained AE loss is given by

$$ $\ell_{\text{cAE, op}}=\ell_{\text{cAE}}+\ell_{\text{cAE, reg}}\quad\text{with} \quad\ell_{\text{cAE, reg}}=\sum_{l=1}^{L}\|\bm{\Phi}^{(l)}_{\text{e}}\|_{ \text{F}}^{2}+\|\bm{\Phi}^{(l)}_{\text{d}}\|_{\text{F}}^{2}+\|\bm{\Phi}^{(l) \intercal}_{\text{e}}\bm{\Phi}^{(l)}_{\text{d}}-\mathbf{I}\|_{\text{F}}^{2}\|( \bm{\Phi}^{(l)\intercal}_{\text{e}}\bm{\Phi}^{(l)}_{\text{d}})^{-1}\|_{\text{F }}^{2},$ (33) $$

with the loss $\ell_{\text{cAE, reg}}$ regularizing the biorthogonality of the AE weights.
Note that this overparametrization does not guarantee the biorthogonality condition, in contrast to the Riemannian approach. In practice, the additional regularization loss $\ell_{\text{cAE, reg}}$ required an additional learning rate, resulting in more possible combinations of hyperparameters.

## Appendix D Training on Acceleration

Fig. 7 illustrates the complete forward dynamics of the proposed RO-LNN trained on acceleration by minimizing the loss (13).

Figure: Figure 7: Flowchart of the forward dynamics of the proposed reduced-order LNN trained on acceleration via (13). The reduction mappings and embeddings of the Lagrangian ROM are depicted in blue and parametrized via a constrained AE with biorthogonal layers. The ROM dynamics are learned via a latent geometric LNN depicted in orange. The mass-inertia matrix is parametrized via the SPD network described in Section 3.
Refer to caption: ROM_network2.pdf

## Appendix E Additional Experimental Details

This section presents additional details on the experimental setup of Section 5.

### E.1 2 2 2 2 - DoF Pendulum of Section 5.1 .

#### E.1.1 Simulation and Data Generation

For this experiment, we use a $2$-DoF pendulum implemented in Mujoco (Todorov et al., ). The pendulum is connected via hinge joints. Both links $i$ are modeled as cylinders of radius $r_{i}=$0.025\text{\,}\mathrm{m}$$, length $l_{i}=$0.4\text{\,}\mathrm{m}$$, and mass $m_{i}=$0.1\text{\,}\mathrm{kg}$$. The simulated environment is entirely free of dissipative effects.
We generate two different datasets $\mathcal{D}_{\text{pend}}=\{\{\bm{q}_{n,k},\dot{\bm{q}}_{n,k},\ddot{\bm{q}}_{n
,k},\bm{\tau}_{n,k}\}_{k=1}^{K}\}_{n=1}^{N}$ consisting of $N$ trajectories with $K$ samples per trajectory. The state evolution is simulated in Mujoco with an RK4-solver with timestep $\Delta t=${10}^{-3}\text{\,}\mathrm{s}$$ over a time interval $\mathcal{I}=[0,T]$\text{\,}\mathrm{s}$$.

Dataset 1: Unactuated pendulum. The initial configurations are randomly sampled from the interval $q_{i}(t=0)\in\left[0,30\right]$\text{\,}\mathrm{\SIUnitSymbolDegree}$$ for each DoF $i=\{1,2\}$. The initial velocities are set to $\dot{\bm{q}}(t=0)=\mathbf{0}$. Each trajectory is recorded for $T=2$\text{\,}\mathrm{s}$$. All models presented in Section 5.1 in the main text are trained on this dataset.

Dataset 2: Sine-Tracking pendulum. Each joint is controlled to follow a reference trajectory via an inverse dynamics torque control law. The reference trajectories are set to follow sinusoidal shapes as
$q_{\text{ref},i}=A_{i}\sin(2\pi f_{i}t+\phi)$, with $\phi_{i}=\arcsin\left(\frac{q_{i,0}}{A_{i}}\right)$, $q_{i,0}\in\left[0,30\right]$\text{\,}\mathrm{\SIUnitSymbolDegree}$$, $A_{i}\in\left[1,30\right]$\text{\,}\mathrm{\SIUnitSymbolDegree}$$, $f_{i}\in\left[\frac{1}{15},1\right]$\text{\,}{\mathrm{s}}^{-1}$$, and $i=\{1,2\}$. $\dot{\bm{q}}_{\text{ref},i}$ and $\ddot{\bm{q}}_{\text{ref},i}$ are obtained as first and second derivatives of $q_{\text{ref},i}$ with respect to time. Initial positions and velocity are set to match the reference trajectories as $q_{i}(t=0)=q_{\text{ref},i}(t=0)$ and $\dot{q}_{i}(t=0)=\dot{q}_{\text{ref},i}(t=0)$
Each trajectory is recorded for $T=$3.5\text{\,}\mathrm{s}$$. Results of models trained on this datasets are presented in App. F.1.

#### E.1.2 Architectures & Training

Geometric LNN architecture.
For all models, the potential energy network $V(\bm{q};\bm{\theta}_{\text{V}})$ consists of $L_{\text{V}}=2$ hidden Euclidean SoftPlus layers of $64$ neurons.
The kinetic energy of all models is parametrized via the mass-inertia matrix of the system.
We consider several architectures for the mass-inertia matrix network $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})$, which differ in their components $g_{\mathbb{R}}$, $g_{\text{Exp}}$, and $g_{\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}$ as follows:

- •
$g_{\mathbb{R}}$: The Euclidean component consists of $L_{\text{T},\mathbb{R}}$ Euclidean layers on the tangent space of the SPD manifold with SoftPlus activations and $64$ neurons per layer.
- •
$g_{\text{Exp}}$: The basepoint $\bm{P}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$ is either set as the origin $\mathbf{I}$, or learned as a network parameter $\bm{P}_{\bm{\theta}}\in\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}$.
- •
$g_{\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}$: The SPD component consists of $L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}$ hidden SPD layers. We consider two different FC SPD layers, namely GyroAI (26), and GyroSpd++ layers (27). Both FC layers can be augmented with a ReEig nonlinearity (29). The regularization constant $\epsilon$ of the nonlinearity is uniformly set as $\epsilon=1e^{-4}$. Note that for the SPD layers, the common notion of number of neurons per layer or layer width is not directly applicable. As the layers act as a transformation $\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}\to\mathcal{S}_{{
\scriptscriptstyle{\mathcal{++}}}}^{n}$ via network parameters $\bm{\theta_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}}$, they can be interpreted as carrying a single transformation per layer, or as transformations of manifold dimension $\text{dim}(\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n})=\frac{n(n+1)}
{2}$.

Baseline LNN architectures.
We compare the proposed geometric LNN against DeLaN, which parametrizes the mass-inertia matrix via a Cholesky network. We consider two variants here: *(1)* The kinetic and potential energy networks are independent as in (Lutter & Peters, ). Both networks consist of $2$ Euclidean SoftPlus layers of $64$ neurons. Hereby, the potential energy network is fully equivalent to the ones employed in the geometric LNN. The kinetic energy network returns a Cholesky decomposition $\bm{L}(\bm{q})$ of the mass-inertia matrix $\bm{M}(\bm{q})$; *(2)* The kinetic and potential energy networks share parameters, i.e., $\bm{\theta}_{\text{T}}\cap\bm{\theta}_{\text{V}}$ as in (Lutter et al., ). The MLP consists of $2$ hidden Euclidean SoftPlus layers of $64$ neurons, and separate output layers return the potential energy and the Cholesky decomposition.

Training on Acceleration Data.
Each training dataset $\mathcal{D}_{\text{train}}$ consist of randomly-sampled data from $30$ trajectories from one of the datasets presented in Section E.1.1. The models are trained by minimizing the LNN loss (7) for $3000$ epochs with architecture-specific learning rates.
If not further specified, both $g_{\mathbb{R}}$ and $g_{\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}$ consist of each 2 layers.

Training on multi-step Integration.
For these experiments, we randomly sample recordings $i$ form $30$ trajectories, which then form the training dataset
$\mathcal{D}_{\text{train}}=\left\{\left\{\bm{q}\right\}_{t_{0,i}}^{t_{0,i}+H
\Delta t},\left\{\dot{\bm{q}}\right\}_{t_{0,i}}^{t_{0,i}+H\Delta t},\left\{\bm
{\tau}\right\}_{t_{0,i}}^{t_{0,i}+H\Delta t}\right\}$. Note that this training dataset depends on the ODE solver’s timestep $\Delta t$ and desired prediction horizon $H$. The models are trained by minimizing the loss

$$ $\ell_{\text{LNN}}=\frac{1}{HN}\sum_{i=1}^{N}\sum_{j=1}^{H}\int_{t_{i}}^{t_{i}+ j\Delta t}\|\dot{\check{\bm{q}}}_{\text{p},i}(t_{i,j})-\dot{\check{\bm{q}}}_{i }(t_{i,j})\|^{2}+\|\tilde{\dot{\bm{q}}}_{\text{p},i}(t_{i,j})-\dot{\bm{q}}_{i} (t_{i,j})\|^{2}+\gamma\;\|\bm{\theta}\|_{2}^{2}$ (34) $$

for $1500$ epochs with architecture-specific learning rates. Notice that the loss (34) is the equivalent of (15) for FOMs. In this category, the featured kinetic energy networks all consist of $2$ hidden layers, and the exponential map layer is implemented with respect to $\mathbf{I}$.

Long-term prediction Experiments. In those experiments, we use both models trained on acceleration data and with multi-step integration. In both scenario, the training datasets contain $|\mathcal{D}_{\text{train}}|=8000$ samples and are trained for $2000$ epochs.

### E.2 Coupled Pendulum ( 16 16 16 16 DoFs ) of Section 5.2.1

#### E.2.1 Simulation and Data Generation

In this experiment, we consider a $16$-DoF coupled pendulum. The pendulum is connected via hinge joints. All links $i$ are modeled as capsules of radius $r_{i}=$0.05\text{\,}\mathrm{m}$$, cylinder length $l_{i}=$0.5\text{\,}\mathrm{m}$$, and mass $m_{i}=$1\text{\,}\mathrm{kg}$$. The simulated environment is entirely free of dissipative effects. We generate a dataset $\mathcal{D}_{\text{pend16}}=\{\{\bm{q}_{n,k},\dot{\bm{q}}_{n,k},\ddot{\bm{q}}_
{n,k},\bm{\tau}_{n,k}\}_{k=1}^{K}\}_{n=1}^{N}$ consisting of $N$ trajectories with $K$ samples per trajectory.
The first $4$-DoF are simulated in Mujoco (Todorov et al., ) with an RK4-solver with timestep $\Delta t=${10}^{-3}\text{\,}\mathrm{s}$$ over a time interval $\mathcal{I}=[0,3]$\text{\,}\mathrm{s}$$, while the last $12$ DoFs are constrained to nonlinear combinations of the first $4$, as described in Table 3.
The initial configurations of the first $4$ DoFs are randomly sampled from the interval $q_{i}(t=0)\in\left[0,30\right]$\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The initial velocities are $\dot{\bm{q}}(t=0)=\mathbf{0}$.

**Table 3: Nonlinear combinations for DoFs 5 to 16.**
| DoF | $f(q_{1},q_{2},q_{3},q_{4})$ |
| --- | --- |
| $q_{5}$ | $q_{3}-\cos(q_{2})$ |
| $q_{6}$ | $q_{1}+0.1\sin(q_{2})$ |
| $q_{7}$ | $q_{4}\cos(q_{2})$ |
| $q_{8}$ | $q_{1}+q_{3}^{2}$ |
| $q_{9}$ | $1.5\sin(q_{2})$ |
| $q_{10}$ | $-q_{4}q_{0}$ |
| $q_{11}$ | $\sin(q_{1})$ |
| $q_{12}$ | $0.4q_{3}q_{4}$ |
| $q_{13}$ | $-0.9q_{1}-q_{2}+q_{3}-2q_{4}^{2}$ |
| $q_{14}$ | $-3\sin(q_{3})$ |
| $q_{15}$ | $-2q_{3}^{2}$ |
| $q_{16}$ | $-0.9q_{1}^{2}$ |

#### E.2.2 Architectures & Training

Our RO-LNN is composed of a constrained AE and a latent geometric LNN. We consider a latent space of dimension $d=n_{0}=4$ to match the known underlying dimension of the high-dimensional pendulum dynamics.
The constrained AE is implemented with $4$ biorthogonal encoder and decoder layers $\rho_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l}}\to\mathbb{R}^{n_{l-1}}$ and $\varphi_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l-1}}\to\mathbb{R}^{n_{l}}$ of sizes $n_{l}=\left[8,16,16,16\right]$. For the latent geometric LNN, both the potential and kinetic energy network consist of $2$ hidden Euclidean layers of $64$ neurons. The kinetic energy network employs an exponential map layer with basepoint $\text{Exp}_{\mathbf{I}}$.

We consider two versions of the RO-LNN. The first model is trained on acceleration data by minimizing loss (13). The second model is trained with multi-step integration by minimizing (15). We consider $H_{\text{train}}=8$ latent integration steps and an integration time constant of ${10}^{-3}\text{\,}\mathrm{s}$. The losses are optimized using Riemannian Adam. For the first model, we use a learning rate of $5e^{-2}$ for the AE parameters $\bm{\Xi}$, $1e^{-5}$ for the LNN parameters $\bm{\theta}$, and a regularization $\gamma=1e^{-6}$. For the second model, the learning rate for the parameters of the LNN is $2e^{-4}$ with $\gamma=2e^{-5}$. Both models are trained until convergence, i.e., for $4000$ and $3000$ epochs, respectively.

We compare our RO-LNN with full-order geometric LNNs. We only report results of the best-performing FOM out of $10$ models. The kinetic energy network consists of $L_{\text{T},\mathbb{R}}=2$ Euclidean hidden SoftPlus layers with $128$ neurons, an exponential map layer with basepoint $\bm{P}=\mathbf{I}$, and no SPD layers, i.e., $L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}=0$. The potential energy network is likewise composed of $2$ Euclidean SoftPlus layers of $128$ neurons.
The full-order LNN was trained for $3000$ epochs. The learning rate for the LNN parameters $\bm{\theta}$ is set to $1e^{-5}$, and $\gamma=1e^{-5}$.

When investigating joint training in Section 5.2, the previously-described RO-LNN trained with multi-step integration is compared to *(1)* a version trained via the overparametrization of biorthogonal weights, and *(2)* a model that was trained separately, by sequentially training the AE and latent LNN. The overparametrized model uses the overparametrized loss (33), with $\ell_{\text{cAE, reg}}$ multplied by a weighting constant $1e^{-5}$.
For the sequentially-trained model, the first $3000$ epochs were trained using only the AE components $\ell_{\text{AE}}$ of the loss (15). For the consecutive $3000$ epochs, only the LNN components $\ell_{\text{LNN}}$ were active of (15). All models were trained on $|\mathcal{D}_{\text{train}}|=24000$ samples, which were randomly sampled from $20$ trajectories. The testing dataset consists of $10$ trajectories.

### E.3 Rope ( 192 192 192 192 DoFs ) of Section 5.2.2

#### E.3.1 Simulation and Data Generation

The rope is implemented in Mujoco via an elastic cable. Over a total length of $1\text{\,}\mathrm{m}$, the rope consists of $i=64$ equally-spaced capsule-shaped masses of $m_{i}=$0.1\text{\,}\mathrm{kg}$$ and $r_{i}=$0.02\text{\,}\mathrm{m}$$. Twisting and bending stiffness are set to $5\text{\times}{10}^{5}\text{\,}\mathrm{Pa}$, and the damping value of the joints is set to $0.1$. As the masses are connected via ball joints, the rope can move in all directions of the workspace.
For training and testing, we generate the datasets $\mathcal{D}_{\text{rope}}=\{\{\bm{q}_{n,k},\dot{\bm{q}}_{n,k},\bm{\tau}_{n,k}
\}_{k=1}^{K}\}_{n=1}^{N}$ consisting of $N$ trajectories with $K$ samples per trajectory. The generalized coordinates of the system are chosen to be the Cartesian positions $q_{i}=\left[x_{i},y_{i},z_{i}\right]$ of each mass’s center of mass in the world frame.
To generate trajectories, we consider a scenario where one end of the rope is fixed to the origin and the other end is controlled to move along a circular trajectory of radius $r$. This mimicks a handheld manipulation task.

For each reference trajectory, a radius $r_{i}$ is sampled randomly from $r_{i}\in\left[0.05,0.4\right]$\text{\,}\mathrm{m}$$. To achieve non-planar motions, the circles are tilted by a random angle $\theta\in\left[-90,90\right]$\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The end of the rope is set to track such a reference at an angular velocity $120\text{\,}\mathrm{\SIUnitSymbolDegree}\text{\,}{\mathrm{s}}^{-1}$ via a PD control.
The state evolution is simulated with an RK4-solver with timestep $\Delta t=${10}^{-3}\text{\,}\mathrm{s}$$ over a time interval $\mathcal{I}=[0,3.8]$\text{\,}\mathrm{s}$$.

#### E.3.2 Architectures & Training

For our experiments on the rope, we consider four RO-LNNs that only differ from one another in the dimensions of the latent space. Each network is composed of a constrained AE and a latent geometric LNN.
The constrained AE is implemented with $4$ biorthogonal encoder and decoder layers $\rho_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l}}\to\mathbb{R}^{n_{l-1}}$ and $\varphi_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l-1}}\to\mathbb{R}^{n_{l}}$ of sizes $n_{l}=\left[64,64,128,192\right]$. We consider latent spaces of dimension $d=\left\{4,6,10,14\right\}$. In all scenarios, the kinetic energy network of the latent geometric LNN consists of $2$ hidden Euclidean layers of $64$ neurons with SoftPlus activations, and an exponential map layer $\text{Exp}_{\mathbf{I}}$. The potential energy network also consists of $2$ hidden Euclidean layers of $64$ neurons with SoftPlus activations.

All models are trained on $H_{\text{train}}=8$ integration steps by minimizing the loss (15) with Riemannian Adam. All models are trained with a learning rate of $8e^{-2}$ for the AE parameters $\bm{\Xi}$, $1e^{-5}$ for the LNN parameters $\bm{\theta}$, with a regularization $\gamma=1e^{-6}$. All models are trained until convergence for $3000$ epochs on $24000$ samples of training data, sampled randomly from $20$ different trajectories. For the testing dataset, we record $10$ full trajectories.

### E.4 Cloth ( 600 600 600 600 DoFs ) of Section 5.2.3

#### E.4.1 Simulation and Data Generation

The deformable thin cloth is modelled in Mujoco as a flexible composite object. $200$ masses of $m_{i}=$0.1\text{\,}\mathrm{kg}$$ are equally spaced over a width of $0.1\text{\,}\mathrm{m}$ and length of $0.2\text{\,}\mathrm{m}$. For this model, the generalized coordinates are given by the Cartesian positions $q_{i}=\left[x_{i},y_{i},z_{i}\right]$ of each masses’ center of mass in the world frame.

For training and testing, we generate datasets $\mathcal{D}_{\text{cloth}}=\{\{\bm{q}_{n,k},\dot{\bm{q}}_{n,k},\bm{\tau}_{n,k}
\}_{k=1}^{K}\}_{n=1}^{N}$ consisting of $N=20$ trajectories during training, and $N=10$ trajectories during testing. The state evolution is simulated with an RK4-solver with timestep $\Delta t=${10}^{-4}\text{\,}\mathrm{s}$$ over a time interval $\mathcal{I}=[0,0.25]$\text{\,}\mathrm{s}$$, resulting in $K=2500$ samples.
The trajectories contain recordings of the cloth falling on a sphere from a height of $0.13\text{\,}\mathrm{m}$ in the center above the origin of the sphere. To vary scenarios, the radius of the sphere is randomly-sampled from $r\in\left[0.02,0.12\right]$\text{\,}\mathrm{m}$$.

#### E.4.2 Architectures & Training

The RO-LNN is composed of a constrained AE with layer sizes $n_{l}=\left[32,64,128,600\right]$ for the $l=1,...,4$ encoder and decoder layers $\rho_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l}}\to\mathbb{R}^{n_{l-1}}$ and $\varphi_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l-1}}\to\mathbb{R}^{n_{l}}$, and a latent space dimension $n_{0}=10$. The kinetic energy network of the latent geometric LNN consists of $2$ hidden Euclidean layers of $64$ neurons with SoftPlus activations, and an exponential map layer $\text{Exp}_{\mathbf{I}}$. The potential energy network also consists of $2$ hidden Euclidean layers of $64$ neurons with SoftPlus activations.

The featured model is trained on $H_{\text{train}}=8$ integration steps via the loss (15) optimized using Riemannian Adam with a learning rate of $5e^{-2}$ for the AE parameters $\bm{\Xi}$, and $2e^{-4}$ for the LNN parameters $\bm{\theta}$, and a regularization $\gamma=2e^{-5}$. THe model is trained for $3000$ epochs.

## Appendix F Additional Experimental Results

### F.1 Learning Lagrangian Dynamics with Geometric LNNs

This section presents additional results on learning the dynamics of a $2$-DoF pendulum with the geometric LNN presented in Section 3.

#### F.1.1 Additional Results on the Experiment of Section 5.1

Trajectory prediction.
Fig. 8 complements Fig. 2-*right* by depicting the long-term position and velocity predictions of the different geometric LNNs considered in Section 5.1. The velocity predictions follow the same trend as the positions: The geometric LNNs outperform the DeLaNs, the geometric LNNs trained via multi-step integration outperform those trained on acceleration data, and the geometric LNNs with SPD layers do not noticeably outperform those employing solely Euclidean and exponential map layers.

Figure: Figure 8: Numerically-integrated position and velocity predictions and ground truth trajectory ( ) for LNNs trained on the unactuated $2$-DoF pendulum dataset.

Energy conservation.
Fig. 9 displays the total energy
of the different LNNs of Fig 8 along the predicted trajectory. As the pendulum is unactuated, the energy levels — that can be learned up to a constant — should be constant over the predicted trajectory. We observe that all LNNs approximately conserve the total predicted energy $\check{\mathcal{E}}$, with the DeLaNs showing the highest fluctuations. Note that, since the predicted trajectories are based on Euler forward integration, we do not expect the total energy of the system to be perfectly preserved.

Figure: Figure 9: Total energy $\check{\mathcal{E}}$ (up to a constant) of different LNNs.

#### F.1.2 Sine-Tracking 2 2 2 2 - DoF Pendulum

We reproduce the experiment of Section 5.1 on a different dataset, namely on the sine-tracking $2$-DoF pendulum dataset introduced in App. E.1.1.

Fig. 2 depicts the acceleration prediction errors for selected architectures.
Similarly as for the unactuated dataset, the LNNs with Cholesky networks are outperformed by the geometric LNNs, and even more so in the low-data regime. In contrast to the unactuated pendulum dataset, the geometric LNNs featuring SPD layers result in slightly lower acceleration errors.

Fig. 11 depict long-term trajectory prediction obtained by integrating the state predictions given by the LNN.
We observe that training the LNNs with multi-step integration is crucial to obtain accurate trajectory prediction for the sine-tracking dataset.
Generally, the results on this second dataset are aligned with those presented in Section 5.1 and validate the effectiveness of considering the geometry of the mass-inertia matrix and of training the models via multi-step integration when learning Lagrangian dynamics.

Figure: Figure 10: Median acceleration prediction error for different LNNs and training set sizes $\mathcal{D}_{\text{train}}$ over $10$ testing trajectories. The models were trained on the sine-tracking $2$-DoF pendulum dataset. Shaded regions represent first and third quartiles.

Figure: Figure 11: Numerically-integrated position and velocity predictions and ground truth trajectory ( ) for LNNs trained on the sine-tracking $2$-DoF pendulum dataset.

#### F.1.3 Influence of SPD layers

Next, we evaluate the influence of SPD layers within the kinetic energy network $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})$ of the geometric LNN.
We consider the network architectures described in App. E.1.2 and train the models with acceleration data. We vary the type of SPD layers, as well as the number of hidden layers $L_{\text{T},\mathbb{R}}$, $L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}$ of the Euclidean and SPD components of the kinetic energy network.
Each model is trained with two random seeds and tested over $10$ unseen testing trajectories from the corresponding dataset. We report the lower prediction error across the two trained models.
The average acceleration prediction errors for the unactuated and sine-tracking datasets are given in Tables 4 and 5, respectively.

Despite the increased performance of geometric LNNs over DeLaNs, the reported prediction errors do not allow us to observe any clear influence of the choice of SPD layers, exponential map basepoints, or Euclidean and SPD layers depths across training dataset sizes. In particular, increasing the depth of SPD layers does not seem to increase the expressivity of the overall SPD network.

We hypothesize that the limited improvements obtained by adding the GyroSpd++, GyroAI, and ReEig layers in our SPD network may be due to practical issues related to the training of these layers. Specifically, the abundance of matrix exponentials and logarithms within the SPD layers tends to cause numerical issues in the optimization procedure. These numerical instabilities also prevent the convergence of several models, as indicated by nan values in Tables 4 and 5. Notice that this mostly occurs for the SPD networks containing ReEig layers.
Potential solutions include implementing analytic layer derivatives to replace unstable autodiff computations, or endowing the SPD manifold with a different Riemannian metric.

Due to the limited performance improvements and increased computational complexity (see App. F.1.5) of the SPD layers), we consider SPD networks $\bm{M}(\bm{q};\bm{\theta}_{\text{T}})=(g_{\text{Exp}}\circ g_{\mathbb{R}})(\bm
{q})$ composed of Euclidean and exponential map components when learning reduced-order Lagrangian dynamics in the experiments of Section 5.2.

**Table 4: Mean and standard deviation of the acceleration prediction errors $\|\ddot{\bm{q}}_{\text{p}}-\ddot{\bm{q}}\|$ for geometric LNNs with different SPD network architectures trained on various training sizes on the unactuated $2$-DoF pendulum dataset.**
|  | $|\mathcal{D}_{\text{train}}|$ |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  | $1000$ | $2000$ | $3000$ | $5000$ | $7000$ |  |  |
| GyroSpd++ | $\text{Exp}_{\bm{I}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $5.95\pm 5.69$ | $5.22\pm 3.78$ | $1.83\pm 2.61$ | $1.86\pm 2.98$ | $\underline{0.92\pm 1.31}$ |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $15.7\pm 19.5$ | $5.15\pm 3.66$ | $4.43\pm 2.96$ | $2.71\pm 3.22$ | $1.52\pm 1.69$ |  |  |
| $\text{Exp}_{\bm{P}_{\bm{\theta}}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $\underline{1.78\pm 2.30}$ | $4.48\pm 2.70$ | $4.45\pm 2.96$ | $1.27\pm 1.47$ | $1.75\pm 2.32$ |  |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $4.37\pm 2.65$ | $5.34\pm 3.72$ | $2.08\pm 3.02$ | $\bm{0.89\pm 1.11}$ | $1.44\pm 2.02$ |  |  |
| GyroSpd+++ReEig | $\text{Exp}_{\bm{I}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $4.66\pm 2.99$ | $1.75\pm 2.33$ | $\underline{1.77\pm 2.40}$ | $1.66\pm 2.58$ | $1.10\pm 1.42$ |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $4.66\pm 2.81$ | $1.72\pm 2.19$ | $1.80\pm 2.52$ | $1.78\pm 2.21$ | $\bm{0.84e\pm 1.25}$ |  |  |
| $\text{Exp}_{\bm{P}_{\bm{\theta}}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $4.60\pm 2.87$ | $4.58\pm 2.93$ | $4.38\pm 2.85$ | $1.43\pm 1.89$ | $1.18\pm 1.93$ |  |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $1.99\pm 2.23$ | $1.78\pm 2.38$ | $1.79\pm 2.52$ | $\underline{1.20\pm 1.77}$ | $1.14\pm 1.31$ |  |  |
| GyroAI+ReEig | $\text{Exp}_{\bm{I}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $2.31\pm 2.33$ | $7.28\pm 8.43$ | $\bm{1.52\pm 2.08}$ | $1.32\pm 1.91$ | $1.13\pm 1.48$ |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $\bm{1.77\pm 2.34}$ | nan | $1.99\pm 2.70$ | nan | $1.30\pm 2.17$ |  |  |
| $\text{Exp}_{\bm{P}_{\bm{\theta}}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $2.76\pm 2.49$ | $\underline{1.71\pm 2.29}$ | $21.2\pm 19.5$ | $1.26\pm 1.66$ | $1.07\pm 1.48$ |  |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | nan | $\bm{1.57\pm 2.01}$ | nan | $1.82\pm 2.48$ | $4.22\pm 2.60$ |  |  |

**Table 5: Mean and standard deviation of the acceleration prediction errors $\|\ddot{\bm{q}}_{\text{p}}-\ddot{\bm{q}}\|$ for different geometric LNNs with different SPD network architectures trained on various training sizes on the sine-tracking $2$-DoF pendulum dataset.**
|  | $|\mathcal{D}_{\text{train}}|$ |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  | 1000 | 2000 | 3000 | 5000 | 7000 |  |  |
| GyroSpd++ | $\text{Exp}_{\bm{I}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $\bm{3.73\pm 4.46}$ | $4.22\pm 4.42$ | $3.95\pm 4.43$ | $3.37\pm 3.54$ | $2.05\pm 2.41$ |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $10.9\pm 8.62$ | $3.95\pm 4.50$ | $3.89\pm 4.49$ | $2.80\pm 2.66$ | $2.18\pm 2.29$ |  |  |
| $\text{Exp}_{\bm{P}_{\bm{\theta}}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $12.0\pm 13.6$ | $\underline{3.58\pm 4.58}$ | $3.70\pm 4.59$ | $3.73\pm 3.95$ | nan |  |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $17.0\pm 21.7$ | $\bm{3.44\pm 4.84}$ | $3.66\pm 4.40$ | $2.75\pm 2.80$ | $\underline{1.49\pm 1.47}$ |  |  |
| GyroSpd+++ReEig | $\text{Exp}_{\bm{I}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $4.59\pm 4.58$ | $4.17\pm 4.43$ | $4.15\pm 4.48$ | $71.0\pm 79.6$ | $5.80\pm 5.23$ |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $13.5\pm 15.3$ | $4.05\pm 4.26$ | $3.95\pm 4.38$ | $5.78\pm 4.91$ | $\bm{1.27\pm 1.25}$ |  |  |
| $\text{Exp}_{\bm{P}_{\bm{\theta}}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $4.29\pm 4.53$ | $4.55\pm 4.47$ | $4.13\pm 4.43$ | $\bm{1.88\pm 1.74}$ | $3.34\pm 3.65$ |  |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $6.13\pm 4.61$ | $3.90\pm 4.42$ | $3.45\pm 4.81$ | $2.95\pm 3.17$ | $2.64\pm 2.82$ |  |  |
| GyroAI+ReEig | $\text{Exp}_{\bm{I}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $6.39\pm 5.73$ | $4.00\pm 4.45$ | $\bm{3.33\pm 4.30}$ | nan | $6.74\pm 5.59$ |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $9.14\pm 9.73$ | $3.79\pm 4.41$ | $3.67\pm 4.60$ | $2.36\pm 2.58$ | $49.1\pm 58.2$ |  |  |
| $\text{Exp}_{\bm{P}_{\bm{\theta}}}$ | $L_{\text{T},\mathbb{R}}=1,L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{<br>\mathcal{++}}}}^{n}}=3$ | $4.02\pm 4.66$ | nan | $4.11\pm 4.34$ | $\underline{2.35\pm 2.08}$ | $1.86\pm 2.13$ |  |
| $L_{\text{T},\mathbb{R}}=L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{<br>++}}}}^{n}}=2$ | $\underline{3.99\pm 4.66}$ | nan | $\underline{3.40\pm 4.55}$ | $2.78\pm 2.68$ | $1.68\pm 1.75$ |  |  |

#### F.1.4 Influence of Multi-step Integration during Training

We assess the influence of the number $H_{\text{train}}$ of steps considered when training our geometric LNNs.
The featured networks are trained as described in E.1.2. For this experiment, the architectures are as described for multi-step integration in section E.1.2.
Table 6 reports the velocity prediction errors averaged over $10$ testing trajectories for models trained on different prediction horizons $H_{\text{train}}$ for each of the $2$-DoF pendulum datasets.
We generally notice a significant decrease in the average error for geometric LNNs trained on longer horizons $H_{\text{train}}$. Ultimately, we select $H_{\text{train}}$ by trading off the increase in performance and the training time. The later is investigated next.

**Table 6: Mean and standard deviation of the $H_{\text{test}}=25$-steps velocity prediction errors for a geometric LNN trained via multi-step integration with different horizons $H_{\text{train}}$ and datasets.**
| Dataset | $H_{\text{train}}=1$ | $H_{\text{train}}=4$ | $H_{\text{train}}=8$ | $H_{\text{train}}=12$ |
| --- | --- | --- | --- | --- |
| Unactuated | $4.61e^{-1}\pm 6.83e^{-1}$ | $6.64e^{-2}\pm 1.73e^{-1}$ | $4.80e^{-2}\pm 1.17e^{-1}$ | $4.80e^{-2}\pm 1.13e^{-1}$ |
| Sine-tracking | $2.32e^{-2}\pm 4.07e^{-2}$ | $8.87e^{-1}\pm 2.05$ | $1.25e^{-2}\pm 1.89e^{-2}$ | $1.36e^{-2}\pm 2.22e^{-2}$ |

#### F.1.5 Runtimes

Table 7 shows the average training times of several LNNs trained on datasets of $|\mathcal{D}_{\text{train}}|=1000$ datapoints for $1000$ epochs. The potential energy network and Euclidean part of the kinetic energy network of all models consist of $2$ hidden Euclidean layers of $64$ neurons. For SPD layers, we only consider GyroSpd+++ReEig layers as they are the computationally most demanding SPD layer combination. The featured model consists of $L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}=2$ of these layers.

Importantly, the well-performing geometric LNN with identity basepoint and $L_{\text{T},\mathcal{S}_{{\scriptscriptstyle{\mathcal{++}}}}^{n}}=0$ trains almost as fast as DeLaN with its Cholesky layers. Instead, integrating GyroSpd+++ReEig layers within the SPD network doubles the training time. Note that similar increases were observed for other combinations of SPD layers.
As expected, training the LNNs with multi-step integration increases the training time.
However, such training leads to increased data-efficiency and improved long-term predictions.

During evaluation, running the well-performing LNN with $2$ Euclidean layers and the $\text{Exp}_{\mathbf{I}}$ layer in combination with numerical Euler integration over 25 steps, takes $0.01191$ secs, which could enable certain online applications.

**Table 7: Training times of various LNNs in seconds.**
| DeLaN | Geometric LNN |  |  |  |
| --- | --- | --- | --- | --- |
| Cholesky | $\text{Exp}_{\bm{I}}$ | $\text{Exp}_{\bm{P}_{\bm{\theta}}}$ | $\text{Exp}_{\bm{P}_{\bm{\theta}}}$, GyroSpd+++ReEig | $\text{Exp}_{\bm{I}}$, $H_{\text{train}}=8$ |
| $88$ | $92$ | $121$ | $189$ | $105$ |

### F.2 Additional Results on the 16 16 16 16 - DoF Coupled Pendulum Experiment of Section 5.2.1

This section presents additional results on learning the high-dimensional dynamics of a $16$-DoF pendulum with geometric RO-LNN.

#### F.2.1 Acceleration vs Multi-step Training

Here, we further evaluate the impact of the losses used to train geometric RO-LNNs. Table 8 shows the average acceleration errors of trajectory predictions of RO-LNNs trained on accelerations or via multi-step integration with horizon $H_{\text{train}}=8$.
We observe that the RO-LNN trained via multi-step integration consistently outperforms the RO-LNN trained on accelerations for all prediction horizons $H_{\text{test}}$. This validates the effectiveness of considering successive states during training to learn the dynamics of physical systems.
Notice that training on acceleration data requires computing the second derivatives in the AE layers, resulting in large-scale matrix multiplications, which not only leads to longer training times, but depending on dimensionality of the system, requires vast amounts of memory. In this sense, training via multi-step integration also avoids the computation of these Hessian matrices.

**Table 8: Mean and standard deviations of the reconstructed and latent prediction errors for different prediction horizons with RO-LNNs trained on acceleration and via multi-step integration.**
|  | $\sum_{j=1}^{H}\|\tilde{\bm{q}}_{\text{p}}(t_{j})-\bm{q}(t_{j})\|^{2}$ | $\ell_{\text{LNN},n}$ | $\sum\ell_{\text{LNN},d}$ |  |
| --- | --- | --- | --- | --- |
| $H_{\text{test}}=8$ | acc | $2.65e^{-1}\pm 3.13e^{-1}$ | $1.95e^{0}\pm 2.14e^{0}$ | $1.70e^{-3}\pm 4.94e^{-2}$ |
| ODE | $\bm{1.86e^{-4}\pm 2.57e^{-4}}$ | $\bm{4.06e^{-2}\pm 8.70e^{-2}}$ | $\bm{9.22^{-4}\pm 3.76e^{-3}}$ |  |
| $H_{\text{test}}=25$ | acc | $2.64e^{-1}\pm 3.12e^{-1}$ | $2.25e^{0}\pm 2.50e^{0}$ | $1.45e^{-2}\pm 4.28e^{-2}$ |
| ODE | $\bm{2.05e^{-4}\pm 2.74e^{-4}}$ | $\bm{1.66e^{-1}\pm 4.85e^{-1}}$ | $\bm{7.78e^{-3}\pm 3.25e^{-2}}$ |  |
| $H_{\text{test}}=50$ | acc | $2.64e^{-1}\pm 3.10e^{-1}$ | $3.26e^{0}\pm 4.62e^{0}$ | $5.31e^{-2}\pm 1.753e^{-1}$ |
| ODE | $\bm{4.03e^{-4}\pm 8.62e^{-4}}$ | $\bm{5.41e^{-1}\pm 1.70e^{0}}$ | $\bm{2.76e^{-2}\pm 1.16e^{-1}}$ |  |

#### F.2.2 Constrained vs Standard Autoencoders

We further investigate the effectiveness of our geometric constrained AE whose parameters are obtained via Riemannian optimization compared to two baselines, namely *(1)* constrained AE optimized via the overparametrization of (Otto et al., ), and *(2)* a regular (unconstrained) AE.
All AEs are trained to reduce a set of position data $\mathcal{D}_{\text{train}}=\left\{\bm{q}_{i}\right\}_{i=1}^{N}$ from $15$ trajectories of the $16$-DoF pendulum to a $d=4$-dimensional latent space.
The constrained AEs are implemented as in Section 5.2.1 and described in App. E.2.2.
The regular AE consists of $4$ Euclidean layers with ReLu activations of sizes $n_{l}=\left[64,64,64,16\right]$ for $l=1,...,4$ encoder and decoder layers $\rho_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l}}\to\mathbb{R}^{n_{l-1}}$ and $\varphi_{\mathcal{Q}}^{(l)}:\mathbb{R}^{n_{l-1}}\to\mathbb{R}^{n_{l}}$.
Our constrained AE is trained using Riemannian Adam with learning rate $5e^{-2}$. The overparametrized and regular AEs are trained using Adam with learning rate $5e^{-2}$ for the overparametrized, and $1e{-2}$ for the regular AE. All models are trained for $2500$ epochs.

Fig. 12 shows the reconstruction errors of the AEs for different sizes. We observe that the constrained AEs consistently outperform the regular one. Moreover, the Riemannian optimization leads to lower reconstruction errors than the overparametrization, especially in low-data regimes.

Figure: Figure 12: Median and quartiles of position reconstruction error for the regular AE ( ), the biorthogonal constrained AE ( ), and the overparametrized constrained AE ( ) for RO-LNN trained on different sizes of training sets.

### F.3 Additional Results on the 192 192 192 192 - DoF Rope Experiment of Section 5.2.2

Fig. 13 shows the AE reconstruction and the RO-LNN predictions ($H_{\text{test}}=25$) for selected DoFs of a test trajectory. We observe that the the RO-LNN accurately models the high-dimensional dynamics of the rope.

Figure: Figure 13: Rope position and velocity predictions ( ) from the RO-LNNs trained via multi-step integration and ground truth ( ). The corresponding AE reconstructions ( ) are depicted for completeness. The model is updated with a new initial condition $(\bm{q}_{0},\dot{\bm{q}}_{0})$ every $0.025\text{\,}\mathrm{s}$ ($H_{\text{test}}=25$).

### F.4 Additional Results on the 600 600 600 600 - DoF Cloth Experiment of Section 5.2.3

**Table 9: Cloth reconstruction and $H_{\text{test}}=25$-steps prediction errors over $10$ testing trajectories for the RO-LNN with $d=10$ trained via multi-step integration with the loss (15).**
| $\|\tilde{\bm{q}}-\bm{q}\|^{2}$ | $\|\tilde{\dot{\bm{q}}}-\dot{\bm{q}}\|^{2}$ | $\ell_{\text{LNN},n}$ | $\ell_{\text{LNN},d}$ |
| --- | --- | --- | --- |
| $2.47e^{-3}\pm 1.15e^{-3}$ | $1.50e^{0}\pm 1.81e^{0}$ | $1.51e^{0}\pm 1.81e^{0}$ | $5.07e^{-4}\pm 1.09e^{-3}$ |

Table 9 shows the reconstruction and prediction errors of our model averaged over $10$ testing trajectories.
Moreover, Figs. 14 and 15 show the AE reconstruction and the predictions of the RO-LNN of selected DoFs of a test trajectory for horizons $H_{\text{test}}=25$ and $H_{\text{test}}=3000$, respectively.
Notice that the ground truth only contains actuation data for the first $0.25\text{\,}\mathrm{s}$ (i.e., $2500$ timesteps). Therefore, we reasonably assume that the cloth remains still for the additional time horizon $t_{\text{aug}}=\left[0.25,0.30\right]$\text{\,}\mathrm{s}$$ and set the ground truth states and torques as equal to their last recorded values.
Fig. 16 displays the predicted latent energy to be compared with the groud-truth energy projected in the latent space. Overall, our results demonstrate the ability of our RO-LNN to infer long-term predictions of complex high-dimensional deformable systems.

Figure: Figure 14: Cloth position and velocity predictions ( ) from the RO-LNNs trained via multi-step integration and ground truth ( ). The corresponding AE reconstructions ( ) are depicted for completeness. The model is updated with a new initial condition $(\bm{q}_{0},\dot{\bm{q}}_{0})$ every $0.025\text{\,}\mathrm{s}$ ($H_{\text{test}}=25$).

Figure: Figure 15: Cloth position and velocity full-horizon predictions ( ) from the RO-LNNs trained via multi-step integration and ground truth ( ). The corresponding AE reconstructions ( ) are depicted for completeness. The dynamics are predicted from a given initial condition $(\bm{q}_{0},\dot{\bm{q}}_{0})$ for $0.3\text{\,}\mathrm{s}$ ($H_{\text{test}}=3000$). The grey-shaded areas indicate the interval during which ground-truth data are extrapolated from the last observation.

Figure: Figure 16: Predicted ( ) and ground truth ( ) latent energy $\check{\mathcal{E}}$ (up to a constant) over time. The grey-shaded areas indicate the interval during which ground-truth data are extrapolated from the last observation.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: Most of this reviewer's questions is dealt with in the weakness section above.

Additionally, how could we connect the satisfaction of the Euler-Lagrange equations of the latent dynamics to that of the original dynamics?
soundness: 3
strengths: The overall structure of RO-LNN is visualized well in Fig. 1, and the authors put great effort into explaining their architectures and learning algorithms that contain lots of equations in the limited space.
summary: This paper investigates the problem of efficiently learning high-dimensional system dynamics from trajectory data. It simultaneously learns a reduced-order model using constrained AE and a dynamics model in the latent space that satisfies Euler-Lagrange equations. Preserving this Lagrangian structure imposes an inductive bias in the latent space to increase the generalization capability of the learned model. The authors proposed a new SPD network architecture to parameterize the mass-inertia matrix in LNN and leveraged Riemannian optimization methods to train the model considering the parameters' geometry.
weaknesses: It seems like the overall idea of learning ROM using constrained AE that preserves the Lagrangian structure is already presented in [Buchfink et al.] cited in the paper. However, this is not clearly pointed out in the paper, and the authors' claims in the introduction and their explanation of the methodology in Section 3&4 somehow mislead the readers to think that the idea is also the main contribution of this paper. This reviewer thinks this paper's contribution on top of the prior work is proposing a new SPD network architecture for parameterizing the mass-inertia matrix and using Riemannian optimization in training to consider the parameters' geometry.

This identification of the contribution should be explained clearly, and the performance of RO-LNN should be compared to the prior works in the experiments.

There are no experimental results that show the benefit of preserving the Lagrangian structure in the latent space. Adding a baseline that learns the model without LNN would help to understand this benefit.

For the proposed architecture of the SPD network in section 3.1, the role of the SPD layers is not explained well, and the experiments show better results without those layers. Justification of having that component is required.

While using Riemannian optimization methods in training RO-LNN is one of the main contributions of this paper, the process is not explained well in Section 3.2.

For eq (13), there is no description of how to compute $\ddot{\check{q}}\_i$. Also, it seems $\tilde{\ddot{q}}\_i = \tilde{\ddot{q}}\_{p,i}$ so that the third term of $\ell_\text{AE}$  should be disregarded, or it should be explained what $\tilde{\ddot{q}}\_i$ is.

In the experiments, the trajectories in fig.2,3 are too cluttered. It might be better to just plot the log-scale difference between the true and predicted trajectories. Also, the plots mislead the readers into thinking that the learned models are used to predict the long-term trajectories while they were fed with ground-truth initial points every 0.025s.

For the cloth experiments, it seems a much shorter timestep size is used for full-horizon prediction (25 steps for 0.025s vs 2500 steps for 0.25s). Then, fig. 6 is not a fair comparison, and this should be pointed out clearly.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. The idea of using Lagrangian principle is usually made to conserve the energy. However, some of the case (e.g., the cloth) does not seem to conserve the energy. How can the authors justify the use of Lagrangian principle in such cases?
2. The error in Table 1 seems to be very high (it reaches $1\times 10^8$). Why can this happen? Normally, a well-setup neural network train should be normalized so the error starts around 1. The error reaches $10^8$ seems to indicate that the error is not normalized. Does this mean the training is not normalized?
3. (Related to weakness #2) Computational complexity in solving the ODE of a system does not only depend on how many dimension it is, but also how stiff the ODE is. Is there any effort in reducing the stiffness of the equation to ensure lower computational demand?
4. How is the roll-out of the NN done? What kind of ODE solver is used? I found in the appendix that the simulation is done using RK4. Is the same method used to roll-out the NN for the ODE loss?
5. Why the number of rolled out time steps $H=8$ quite small in computing the ODE loss? What is preventing it from using higher number of time steps?
6. How was the time step $\Delta t=10^{-3}$ chosen? How would choosing larger time step affect the roll-out prediction error?
soundness: 2
strengths: 1. The paper presents experimental results on increasing complexity of the case
2. It shows better results compared to Lagrangian Neural Network
summary: The paper presents an idea of a neural network to learn the dynamics of a system using Reduced Order Model (ROM). ROM is achieved by mapping the state-space dimensions into a lower dimension and learn the dynamics in the reduced lower dimensional space. The authors proposes a neural network that learns the mass matrix and the potential function in the lower dimensional, and compute the acceleration based on Lagrangian principle. The model is trained in 2 ways: (1) learning the acceleration, and (2) learning the trajectory (and velocity) in a roll-out dynamics. The results show significantly better results than Lagrangian Neural Networks.
weaknesses: Major:
1. It is unclear where the novelty of the idea lies and what advantage it brings with respect to the previous works. The use of auto-encoder to reduce the dimensionality of the system is not a new idea. In fact, several papers in the same field like Hamiltonian Neural Network, Lagrangian Neural Network also used auto-encoder to learn dynamics of systems from images. If I understand it correctly, the novelty of the idea lies in a neural network learning both the mass matrix and potential for the Lagrange term. However, it is unclear what benefits does this new formulation of Lagrange term has in comparison with, for example, DeLan and LNN.
2. Reducing the dynamics to lower order dimension does not mean lower computational demand. It is possible, for example, that the lower dimensional dynamics has stronger stiffness than the original dimension. The paper does not present any results comparing the wall-clock time of solving ODE in the original dimension vs solving learned ODE in the lower dimension using the same hardware (e.g., CPU and GPU) to be convincing that their method really achieve computational reduction.
3. Lack of comparisons with other architectures. It seems like the authors only comparing the experiment with LNN. There are a lot of other works in this area that the authors can compare their work with. For example, Hamiltonian Neural Networks (https://arxiv.org/abs/1906.01563), dissipative HNN (https://arxiv.org/abs/2201.10085), dissipative SymOden (https://arxiv.org/abs/2002.08860), Constant of motion network (https://arxiv.org/abs/2208.10387), etc.

Minor comments:
1. Please avoid using the unscientific notation like $1e^{-5}$. Instead, please use notations like $1\times 10^{-5}$.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 3
questions: see the weakness part.
soundness: 3
strengths: The strengths are listed as followed
- Utilizes geometry and physics as inductive biases, leading to a physically consistent model to improve the interpretability. 
- Effectively learns reduced-order dynamics for high-dimensional rigid and deformable systems with accurate long-term predictions.
- Joint optimization of latent space and reduced dynamic parameters enhances performance, particularly in low-data regimes.
- Incorporates intrinsic geometry of the problem, aiding in the preservation of energy and interpretability of the model.
- Outperforms standard autoencoders by employing a constrained auto-encoder trained with Riemannian optimization.
summary: The paper proposes a geometric network architecture that learns reduced-order Lagrangian dynamics for high-dimensional systems, using a Riemannian approach to jointly learn a structure-preserving latent space and the associated low-dimensional dynamics. This method efficiently predicts long-term dynamics for rigid and deformable systems while improving interpretability based on physics.
weaknesses: ### Contribution and Originality

- I question the novelty of this paper. The proposed approach closely resembles existing methods described in [1, 2]. Additionally, the authors have not cited reference [1], which appears to be highly relevant to their work. 

- The contribution of this paper is overstated, and the true contribution appears to be quite incremental. Based on my understanding, the theoretical foundation relies heavily on the published work in [2], with many of the equations directly derived from [2]. There seem to be no novel contributions to the methods or framework, as the use of an autoencoder for reduced-order modeling was already proposed in [2]. This paper primarily appears to be an implementation of the original idea from [2].

- For the learning framework, this paper replaces the Galerkin method (subspace method) with an autoencoder structure for learning Lagrangian dynamics. However, the autoencoder representation is fundamentally similar to the Galerkin (subspace) method in its underlying principles. The primary distinction from reference [1] lies in the representation of the mass-inertia matrix as a symmetric positive definite (SPD) manifold within the configuration space, whereas the loss functions and implementation details remain largely consistent with those in [1]. The authors should clearly elucidate the key differences from prior work and provide a justification for why reference [1] was not cited.


### Comparison

- There is no comparison with other existing algorithms [1, 3]. The authors claim scability of its algorithm, but the current work only tests up to 600 DoFs using deep learning. In contrast, the existing work in [1] conducted experiments with 259,744 DoFs, even using Matlab. This raises questions about the claimed scability and performance improvements.

### Computational Complexity

- The computational complexity is significantly higher compared to [1], mainly due to the learning of the SPD manifold. It does not seem necessary to optimize such a complex manifold for high DoF cases, as the constant symmetric positive-definite surrogate used in [1] performs adequately in a complex scenario. 


[1]. Sharma, Harsh, and Boris Kramer. "Preserving Lagrangian structure in data-driven reduced-order modeling of large-scale dynamical systems." Physica D: Nonlinear Phenomena 462 (2024): 134128.

[2] Buchfink, Patrick, et al. "Model reduction on manifolds: A differential geometric framework." Physica D: Nonlinear Phenomena 468 (2024): 134299.

[3] Harsh Sharma, Hongliang Mu, Patrick Buchfink, Rudy Geelen, Silke Glas, and Boris Kramer. Sym- plectic model reduction of Hamiltonian systems using data-driven quadratic manifolds. Computer Methods in Applied Mechanics and Engineering, 417:116402, 2023. doi: 10.1016/j.cma.2023. 116402.

---

confidence: 2
contribution: 2
initial_rating: 6
presentation: 2
questions: Q1. Please refine the notation system in the paper. You may remove some unnecessary notations and descriptions. See my Weakness points 3, 4, and 7.

Q2. The background and design motivation should be emphasized for readers unfamiliar with Lagrangian systems. See my Weakness points 1, 2.

Q3. The writing logic should be improved. See my Weakness points 5, 6.

Q4. The computational cost should be illustrated and reported. See my Weakness point 8. 

Q5. The result presentation can be improved. See my Weakness points 9 and 10.
soundness: 3
strengths: The author proposes a dimensionality reduction framework with physical consistency and geometry preservation. The writing is clear and the experiment demonstrates the effectiveness of the model.
summary: The paper proposes a reduced-order learning model with the preservation of the state geometry for high-dimensional Lagrangian systems.
weaknesses: 1.	In the Introduction, for readers unfamiliar with the high-dimensional Lagrangian system, the author should provide some references and emphasize their importance. The author is also suggested to provide some quantitative descriptions, e.g., the range of the dimensionality.

2.	Is preserving the geometry important between FOM and ROM? Can the author summarize what the benefits are of preserving the structure in physical systems and control? 

3.	What’s the difference between Eq. (3) and Eq. (4)? For example, in Eq. (3), you state $X|\_{\gamma(t)}\in \mathcal{T}\_{\gamma(t)}\mathcal{T}\mathcal{Q} $. In Eq. (4), however, the set is $\mathcal{T}_{\gamma(t)}\mathcal{Q}$. Why do smooth trajectory lie in $\mathcal{Q}$ rather than $\mathcal{T}\mathcal{Q}$, where is $\dot{q}$?

4.	What is $V(q;\theta_V)$ in Line 162? Since there are massive notations, the author is advised to complete the proofreading. 

5.	It’s unclear of the transition between Section 2.3 and 2.4. For example, to fix what problems in Section 2.3, should we introduce LNN in Section 2.4?

6.	The advantage of the proposed SPD network should be better illustrated. The authors claim that their model makes use of Riemannian geometry. However, it seems that other methods can also maintain this constraint. For example, Cholesky decomposition can also make sure the output belongs to SPD space. 

7.	Notations in Eq. (8) are unclear to the reader. What are $l$, $h$, $x$, and PT? Can the author align them with the previous notations? 

8.	It seems weird to state that computing the second-order derivative in Eq. (13) is expensive, but computing the integration of multi-step in Eq. (14) is cheap.  

9.	The legends in Fig. 2 are hard to align with the curves. Please refine this figure. 

10.	Why is the error in Table 1 very high? I am not sure if this is a validated result.

Meta-review: **Summary** This paper provides a geometric approach to learning high dimensional dynamical systems.  The method involves an order reduction autoencoder model which maps system states to a lower dimensional reduced manifold. The order reduction network has layers mathematically constrained to ensure the reduced manifold is embedded and is optimized using Reimannian gradient descent. Dynamics are learned in the lower dimensional latent space with a  Lagrangian neural network approach in which an SPD network parameterizes the mass-inertia matrix and an MLP parameterizes the potential function. Together they define equations of motion in the reduced manifold.         
 
**Strengths** This paper presents an effective method leveraging physical and geometric inductive biases for improving physical fidelity in dynamics modeling and obtaining high accuracy long-term predictions even in systems with many degrees of freedom and unknown dynamics.  The paper is clear.  The empirical evaluation shows the method improves upon previous methods with similar goals such as Lagrangian Neural Networks and L-OpInf (Sharma & Kramer).

**Weaknesses** The main shared concern of reviewers was that the work showed limited novelty relative to previous work, especially Buchfink et al and Sharma & Kramer.  The authors expanded extensively in the rebuttal and revision on the novelty and advantages relative to these and other works, convincing most reviewers.  This leads me to believe that this issue was partly a matter of the writing and evaluations in the initial paper which did not focus enough on clarifying and empirically justifying these differences relative to previous work.  That said, while most reviewers increased their scores to accept, the consensus belief is still that the significance of the contribution is not that larger. Additional concerns which were largely addressed during the rebuttal include the need for more baseline comparisons and a request to show the time complexity.

**Conclusion** The authors largely addressed most reviewer concerns and showed a reasonable contribution relative to previous work.

Decision: Accept (Poster)