A Riemannian Framework for Learning Reduced-order Lagrangian Dynamics
Authors: Katharina Friedl
Noémie Jaquier
Jens Lundell
Tamim Asfour
Danica Kragic
Venue: ICLR 2025 Poster (2025)
Keywords: physics-inspired networks
dynamics learning
model-order reduction
Riemannian geometry
deformable objects
Primary area: learning on time series and dynamical systems
TL;DR: a physics-inspired geometric model to learn dynamics of high-dimensional systems via structure-preserving model-order reduction

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 non-linear 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.

Published as a conference paper at ICLR 2025
A RIEMANNIAN FRAMEWORK FOR LEARNING REDUCED-ORDER LAGRANGIAN DYNAMICS
Katharina Friedl' Noémie Jaquier' :2 Jens Lundell'
Tamim Asfour² Danica Kragic'
1 Division of Robotics, Perception, and Learning 2 Institute for Anthropomatics and Robotics
KTH Royal Institute of Technology (kEriedl, jaquier, jelundel, dani fanijakth.sse seth. asfourkkit.
Karlsruhe Institute of Technology (KIT)
ABSTRACT
By incorporating physical consistency inductive bias, deep neural networks display increased generalization capabilities and data efficiency in learning non- linear dynamic models. However, the complexity of these models generally in- creases with the system dimensionality, requiring larger datasets, more complex deep networks, and significant computational effort. We propose a novel geo- metric network architecture to learn peysscalll-consistern reduced-order dynamic parameters that accurately describe the original high-dimensional system behav- ior. This is achieved by building on recent advances in model-order reduction and by adopting a Riemannian perspective to jointly learn a non-linear structure- preserving latent space and the associated low-dimensional dynamics. Our ap- proach enables accurate long-term predictions of the high-dimensions dynamics of rigid and deformable systems with increased data efficiency by inferring inter- reduced I Lagrangian models.
pretable and
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 net- works (Lutter et al. 2019 Greydanus et al.] [2019] | CCanmer et al., 2020 Lutter & Peters (023) and Gaussian processes (Tnaaa et al.) 2022 Evangelisti & Hirche, 2022). Several physically-consister
models were extended to capture non-conservative forces and contact mechanics (Hochlehnert et al.)
2021 Zhong et al. 2021). and to ensure longer- tergerrttermaeee stability of the identified dynamics by including
differential solvers (Chen et al. 2018) into the training loop (Finzi et al.] 2020; Zhong & Lconard
2020). Model performances were also improved by representing the dynamics using coordinates
adapted to the physical system at hand (Finzi et al.) 2020 Celledoni et al. 2023) Duong & Atanasov
2021]. Altogether, gray-box approaches are physically consistent, versatile, more data-efficient, and display increased generalization capabilities compared to black-box methods (Lutter & Peters) [023]] Greydanus al. 2019). Still, they have, so far, only learned the dynamics of low-dimensional sys- tems, i.e., typically 2-5 dimensions. Learning the dynamics of high-dimg sigh-dimensional systems, such as fluid flows, continuum mechanics, and robots, arguably remains an open problem. This is due to their increasing complexity, requiring more complex network architectures and more training data.
Predicting trajectories of high-dimensional systems is also notoriously difficult due to the com- putational 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 reduced-order model (ROM) of given high-dimensional dynamical system or full-order model (FOM) with known dynamics (Schilders et al. 2008). MOR is typically achieved by projecting the FOM onto a lower-dimensional space via linear (hheffry
et al. 2019; Farhat et al. 2015) Carlberg et al. 2013) or nonlinear (Sharma et al. 2023 Barnett & Farhat 2022) mappings. In particular, Autoencoders (AEs) were shown be well suited to

---

Published as a conference paper at ICLR 2025
extract descriptive nonlinear reduced representations from data (Lee & Carlberg 2020: Buchfink
iet al. 2023 Otto et al. 2023). Recent approaches investigated the preservation of physics-induced
geometric structures, 1.c. Lagrangian or Hamiltonian structures, of FOMs within projection-based
ROMs (Carlberg et al., 2015; Lee & Carlberg 2020 Otto et al. 2023 Hesthaven et al. 2022). This structure awareness enables the derivation of shared theoretical properties such as energy conser- vation and stability preservation (Buchfink et al. 2023), as well as the design of low-dimensional control strategies on the ROM (Lepri et al.) 2024). Buchfink et al. (2024) recently unified vari- ous structure preeeerviggee MOR techniques by adopting differential differential-geomeere -eoomtrric perspective on the problem. In this framework, ROMs are defined on embedded submanifolds of the high high-dimensional Riemannian or symplectic manifold associated with a given Lagrangian or Hamiltonian system.
Despite their advantages, MOR methods are intrusive in that they assume known high-dimensional dynamic parameters which are generally difficult to obtain. [Sharma & Kramer (2024) presented a novel non-intrusive approach that identifies low-dimensional dynamic parameters structure- preserving linear subspace obtained from high- high-dimensional state observations. However, this ap- proach is limited to systems that are linearly reducible and displays limited expressivity by con- straining the range of the ROM parameters. In this paper, propose more general and expressive approach in the form of a physics -nnsiiiee geometric deep network that learns the continuous time dynamics of high-dimensional systems. Our first contribution is adopt differential geome- try perspective and jointly learn (1) a structure-preserving non-linear reduced representation of the high-dimensional generalized coordinates and (2) the associated low-dimensional dynamic param- eters. Specifically, our approach leverages a constrained AE (Otto et al. (2023) to learn a latent embedded submanifold in which physically-consistenn low-dimensional dynamics are learned with a Lagrangian neural network (LNN) (Cranmer et al., 2020 Lutter & Peters 2023). Our approach
differs from (Greydanus et al. 2019 Zhong & Leonard [220 Botev et al.] (021) 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, intro- duce 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 duced 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, rope, and a thin cloth. Our results demonstrate that it efficiently learns reduced-order dynamics leading to accurate long-term predictions of -igh-dimennional systems.
models as both the ROM and the FOM conserve the re-
BACKGROUND
This section introduces the mathematical preliminaries and network architectures for learning re- duced Lagrangian dynamics. We adopt a cetrrential-gemmettie perspective to simultaneously inves- tigate the use of nonlinear MOR projections and the structure preservation of Lagrangian systems.
2.1 LAGRANGIAN DYNAMICS ON THE CONFIGURATION MANIFOLDD
We consider an n-degrees-of -freeedooo (DoF) mechanical system whose configuration space is iden- tified with an n n-dimensional smooth manifold Q with a simple global chart. Velocities q at cach configuration
Q lie in the tangent space Ta Q, n-dimensional vector space composed of all vectors tangent to Q at q. The disjoint union of all tangent spaces T Q forms the tangent bundle TQ, smooth 2n-dimensional manifold. The configuration manifold can be equipped with Riemannian metric, i.e., umoothly-varyee which, given choice of local coordinates, equals the system's 's mass-inertia matrix M (q).
inner product acting on T Q. We consider the kinetic-energy metric,
A I Lagrangian system is a tuple (2, L) of Riemannian configuration manifold Q and a smooth time- independent Lagrangian function L: TQ R. The Lagrangian function is given by the difference between the system's 's kinetic [qq, q) and potential (q)) energies as L(q. q) = l'qq. q) - V(g) = g (q) V(q). Following the principle of least action, the equations of motion are given by the Luler-Lagrange equations d (aC) ac
oq T with generalized forces T. These equations are
de
d M(g) + c(g. 9) + + r, with c(q. q) M(g))
1 (q (QTM(a)d))T, (1)
dt
2 2g
2

---

Published as a conference paper at ICLR 2025
and g(q) oVia) öq where c(q. Q) represents the influence of Coriolis forces. Given time interval
I [tv. ty], trajectories 7 I - TQ: + q(t) of the Lagrangian system obtained by solving the initial value problem (IVP)
d >/. = X 7(t) dr
4(0) M ((0)(r-c(a.)) -(a.0))-a00) -(a.))))aa))))
E T
(2)
((to)= = (
E ETQ,
with X (y(t) X the vector field defined by the Euler-Lagrange equations. Note that, for Lagrangian systems, trajectories are lifted curves on the tangent bundle TQ
2.2 MODEL ORDER REDUCTION
MOR methods consider a high-dimensional dynamical system, i.e., a FOM, with known dynamic parameters, which generates smooth trajectories 7 IZ M characterized by an IVP of the form
d
d de
x x M, tEI,
(3)
(tt)
= Yo E M.
where X (x(e) is smooth vector field defining the evolution of the system, so that 7(1) M. Note that we have M = TQ in ( The goal of MOR is to accurately and efficiently approximate the set of solutions S (t()) E MM S M of the IVP G To this end, MOR methods learn
reduced manifold M with dim dd < dim(M)
and ROM $ $7
X 4c0) E Tra
Following Buchtink et al. (2024), we identify M with embedded submanifold ((M) S M via smooth embedding 3 MM + M. The reduced initial value 500 p(7o) and vector field 4(6) = dp r(t) X ((1)) are defined via the point and tangent reduction maps M + M and
X
dplæ M - M associated with yo, where a .MM The reduction maps must satisfy the projection properties
(4) Trajectories of the original system are then obtained via the approximation r(t) ((+(t)). In this paper, we consider the case where the high-dimensional dynamics are unknown. We lever- age structure-prexerving MOR methods to learn a non-linear reduced representation Q of the high- dimensional Lagrangian configuration space Q along with Lagrangian ROM (Q. C).
pop=id and dpla(a) delæ = idTT.M Vz M.
A key aspect of MOR is the construction of the embedding &o and corresponding point reduction p. AEs well suited for representing these nonlinear mappings due their expressivity and scalabil- ity (Lee & Carlberg 2020 Otto et al. 2023). An AE consists of an encoder network p R -> Rdd that maps a high-dimensional vector æ R" latent representation a decoder network : IRd + R" that reconstructs an approximation of the original data. Given dataset (æ; (x=);, the AE parameters v are trained by minimizing the reconstruction error
g with < n, and
N
AE V
Σ N
1
Il p(a.) - x; 12:
(5)
Notice that the projection properties hold approximately when the loss (5 is small. However, in this paper, we leverage constrained AE (Otto et al.) 2023) that guarantees which is essential when learning reduced representations of high-dimensional Lagrangian systems.
2.3 LAGRANGIAN NEURAL NETWORKS
In contrast to MOR, LNN consider low oiwddmeesioonal systems with unknown dynamics that we aim to learn. 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-plause
this context, DeLaN (Lutter et al. 2019 Lutter & Peters 2023) proposed model the kinetic energy (via the mass-inertas matrix) and the potential energy of a Lagrangian function L as two networks M (q: 07) and V(g: 0v) with parameters 0 (0r, 0v). Trajectories of the learned dynamical
trajectories.
3

---

Published as a conference paper at ICLR 2025
system then obtained via the equations of motion ① by solving the IVP 2 Given set of N observations (q:: 9i Ii: Ti i= the networks are trained to minimize the loss
N
NN (0)
Σ N
1
If(q. q., T: @: 1l
(6)
1=1
where f(q.q, M- (q; 0T)(T c(g, g(q: 0v) and X is regularization constant.
Following this approach, the total energy E T + V(g) = M q + V(q), of un- forced systems is conserved. Next, we propose an enhanced geometric version of DeLaN, which we then use in Section to learn the dynamics of a Lagrangian ROM.
3 GEOMETRIC LAGRANGIAN NEURAL NETWORKS
The mass-inertia matrix M (q) acts as the Riemannian metric on the configuration manifold Q for Lagrangian systems with quadratic kinetic energy structure. As such, (()) belongs the manifold of symmetric positive-definite (SPD) matrices S (hhatia 2007 Pennec et al. 2006). Existing LNÑS enforce the symmetric rusitve-definiteness by predicting the Cholesky decomposition L(q) with a Euclidean network and describing the mass-inertia matrix as M = LL However, this parametrization overlooks the Riemannian geometry of both Cholesky and SPD spaces (see also App. B). 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 q + (q): Q - S" where S = f Sym" Σ > 0} the Riemannian manifold of SPD matrices and Sym" is the space of symmetric matrices (see App. A.1 for 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 the SPD manifold by leveraging its gyrovectorspace structure (López et al. 2021 Nguyen 2022 Nguyen al. 2024) However, these networks consider that both inputs and outputs belong S
We propose a novel SPD network (q; OT) = (gsg, gExp o (R)(g) with three consecutive compo- nents, as illustrated in Fig. [Heft. First, the input configuration q E Q fed to a standard Euclidean multilayer perceptron (MLP) IR : F + m=a++))/2 whose output elements are identified with those of a symmetric matrix U Sym". Second, U is interpreted as an element of the tangent space TpS" ++ and mapped onto the SPD manifold via an exponential map layer X Exp (U)) (U). Note that this type of layers is commonly used in hyperbolic neural networks (Shimizu al. 202
Third, X is fed to a of Lypa layers IS" S + S",,, which outputs the mass-inerta matrix M E ST We consider and compare two different fully-connected (FC) SPD layers which are analogous the FC layers in MLPs -, namely the gyrocalculus-basee (GyroAI) layers of (López et al. 2021 Nguyen| 2022 and the gyrospace hyperplane-based (GyroSpd._.) layers of (Nguyen 2024). Both can be integrated with the ReFig layers (Huang Gool 2017) which are the SPD counterparts to ReLU-activations (see App. C for details on the SPD layers). In total, the pa- rameters Or of our SPD network consist of the MLP parameters ÔT,R, and the SPD layer parameters 0r.s%, Sn X ST including the basepoint P of the exponential map.
al.
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 M ((; 0r) and the potential energy via a standard Euclidean MLP (iq: 0y). The network parameters 0 = {ØT. 0v) } are obtained for classical LNNs by minimizing 6 However, as several parameters of the SPD network belong to S" their geometry must be taken into account during training. To do so, we leverage Riemannian optimization (Absil et al.) 2007 Boumal 2023) to solve a problem of the form miDæEM ((x) on a Riemannian manifold M. where f' is the loss to minimize and æ E M is the optimization variable. Conceptually, each iteration step in a first-order (stochastic) Riemannian optimization method consists of the three successive operations
(7) First, the update n, E M is computed as a function of the Riemannian gradient grad E at the estimate æ; and of 7. transpoere previous update 7)-1 E - M to the new Ta M
ni + (grad ((I)) Øt+1 + Expa, (-ar7)) T: +- PT (Tt).
4

---

Published as a conference paper at ICLR 2025
Geoenctri LNN
condition
pretictions
Pa
8ym
over time !
(4)
R..
ges
4: 4g
PC dpcla
92 dpcl
f( 9.
4p
X(g)
Ufa) Repp
dpoli
(... 91)
Lagrangian ROM
cAR
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 strained 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).
(see App. for the relevant Riemannian operations). Then, an update of the estimate I; is obtained via projecting the update Th, that is scaled by a learning rate learning rate ag onto the manifold the tangent space
with the exponential map. Finally, the current update is traaslee-trnnsported of the updated estimate to prepare for the next iteration. Note that the function h is determined by the optimization method. We the Riemannian Adam (Becigneul & Ganea 2019) implemented
in Geoopt (Kochurov et al. (0200) to optimize the geometric LNN parameters. Notice that M is defined as a product Euclidean and SPD manifolds optimize the parameters 0.
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 power-dimensionss surrogate dynamic model. Building this assumption, propose to reduce the problem dimen- sionality by learning a latent dynamical system that preserves the original Lagrangian structure. The latent space and the reduced dynamics parameters are inferred jointly via struccuue-preservving 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. with high number of DoFs with observations {q:, Ôt, qi, T ) and not nigh-ddmensional represen- tations of inherently oow-dimesioona state spaces as would be the case with images.
Notice that we consider physical systems
4.1 LAGRANGIAN REDUCED-ORDEL MODEL
Preserving the properties of the original FOM is crucial ensure that the learned ROM displays similar behaviors. Therefore, we leverage the geometric framework introduced in Buchfink et al. 2024) and learn a reduced Lagrangian (Ộ, L) via structure- structure-preservinne preserving MOR. As Lagrangian trajec- tories lifted curves (t) (q(t), q(t)) the manifold be reduced is the tangent bundle TQ We define the lifted embedding ş TQ - TQ for smooth embedding se Q + Q as the pair p(q. q) = (Pp(4), dppla4))
(8)
with depll T - Tp)(g) 2 denoting the pushforward, or differential, of yo at q. Analogously, we define a reduction map p(q, q) as the pair (po(g), d of point and tangent reductions
Po: 2 Q and dpelg Ta 2
Given yp and the observed states {q:: q: 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
Ł= L s.
The Fuler-Lagrange equations of the reduced Lagrangian yield
M M( + č(g. à) + g()) =
(10)
In the intrusive case (Buchfink et al. 2024), the reduced parameters are given as function of the known high-dimensional dynamics M M(g) = de M(q) depla- ğ(q) ddell g(q).
5

---

Published as a conference paper at ICLR 2025
=ddll т, and č(q. q) is computed function of (q) as in I Instead, assume known dynamics and propose learn the reduced parameter with a geometric LNN Specifically, parametrize the reduced mass-inertia matrix and the reduced potential energy of E as two net- works (q: 0f) and (g: 0²) with parameters 0 [0,, 0v). The ROM (TO) can then be used to efficiently compute reduced trajectories (tt) as solutions 7(t) se (((((t) of the FOM ( It is worth emphasizing that the preservation of the Lagrangian structure in the ROM guarantees the conservation of the reduced total energy & along the solutions ((t) of (IO) and the corresponding
image curves ((((t)) as & = E o sp This property allows us to learn physically-consiste reduced dynamics. Next, discuss how learn fo and PQ-
4.2 LEARNING THE EMBEDDING AND POINT REDUCTION
For increased expressivity, we parametrize the point reduction pQ and embedding Pe of our reduced Lagrangian as the encoder and decoder of an AE. To ensure that the reduction map p(g. q) satis- fies the projection properties 44). leverage the constrained AE architecture introduced by Otto et al. (2023). The encoder and decoder networks are given as a composition of feedforward layers 2
(1)
(L) and ip = o. with Po
(1) R - R and (1) R R 4
Pe P2 .... Po
d = no .... S n n, where d denotes the latent space dimension. Notice that the pushforward depla a and tangent reduction deola are the differentials of the encoder and decoder networks.
The key to guarantee the projection properties (4) lies in the specific construction of layer pairs as po (q))) = (vT(qa) bi)) and
(1)
11) (q
)) = 2i0 + (q"t-1)) + bb,
(11)
where (i), 1) and (a-.a-) are pairs of weight matrices and smooth activation functions, respec- tively, and b, are bias vectors. By constraining the pairs to satisfy / o, Ia and a_ 0,
each layer satisfies (1) (1) idg and the constrained AE fullfills (4) Otto et al. 20033 adhere to the first constraint by defining pairs of biorthogonal matrices (Φ), 44) via an over- parametrization resulting additional loss terms. Instead, we adopt geometric approach that accounts for the Riemannian geometry of biorthogonal matrices. Specifically, we minimize the AE reconstruction error (5) via Riemannian optimization (see Equation (7)) by considering each pair (L,, 1) as elements of the biorthogonal manifold Bmt, = { V/ R x R"
Pe % po
T I (see App. A.2 for background on the biorthogonal manifolds and associated op- erations). 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, invert- ible activation functions defined (Otto et al. 2023 Equation 12). Additional details on the AE architecture, including activation functions, and layer derivatives, are provided in App.
4.3 MODEL TRAINING
Finally, we propose to jointly learn the parameters E 1 Vi, bi of the AE and 0 (0 0r,sq_ (0) of the latent geometric LNN. We consider two losses, both of which minimize the AE reconstruction error, the latent LNN loss, and joint error on the reconstructed predictions.
Training on acceleration. Given a training set {gi; qi, Ts (=1" the acceleration loss is
N
PROM.
4 4:11
1|: d:1l" 1/
iqp
+X10013.
(12)
N
EAA
fINN,
T.. d: ÿ: are the reconstructed position, velocity, and acceleration obtained by successively applying the lifted point reduction and embedding, and its derivatives, and q is the reduced acceleration. The latent acceleration predictions : qp.i obtained from the reduced equations of motion (00) as
To = f( á + 0) = - 0+) č(g 8) g(:: 0c)),
(13)
The acceleration of the FOM is then computed as 7. de + dl pel à shown in Fig
Training with multi-step integration. The loss (22 requires computing the second derivatives of the embedding and point reduction, i.e., the Hessian matrices of the AF, which result in significant

---

Published as a conference paper at ICLR 2025
computational efforts increasing with the dimension of the FOM. 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 (13) before decoding. This is achieved via H Euler forward steps with constant integration time At. Given sets of observations {q, (I,), q: (I,), (I,)
over intervals Ij = [t., t. + H Htt,, the multi-step loss is
PROMODEE
ΣΣ HN
$ (tr.s) - q (fr (t()) q. (ta.)))
(14)
+ lggpe(tes) - (c(())))) (tps(tes) (((tn)) +y 1011.
with latent velocity predictions ap.:(t..)) Pai f (t:(t), (t), + (t); 0), velocity reconstructions
( = doo (t.(t;))) and tij t, jAt. Fig. iillustrates 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 learn the dynamics of three simulated nigh-dimensional systems: 16- DoF pendulum, a 192-DoF rope, and a 600-DoF thin cloth. Our experiments demonstrate the ability of our approach learn reduced-order dynamics, resulting accurate long-term predictions. More- over, 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. Eop.Ffoo for each experiment. Additional results are provided in App. G A video and source code are available at htttttss//siteee google. com www.ooovwewwwovvoco
lagrangiann
5.1 LEARNING LAGRANGIAN DYNAMICS WITH GEOMETRIC LNNS
We start by evaluating the long-term prediction accuracy the geometric LNN introduced Sec- tion Jon 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 Cholesky network. We consider two variants, where the kinetic and potential energy networks are indepen- dent (Lutter & Peters 2023) or share parameters (0r n Ov) (Lutter al. 2019). Moreover, evaluate the performance of different architectures for our SPD network M (q; 0T). First, we com- pare exponential map layers, where the basepoint P S is either set as I or learned parameter. Second, we evaluate the influence of the SPD layers by comparing performances without any SPD layers, and with GyroAI+ReEig, GyroSpd +. CcyroSpd,,+ReEi. layers.
Fig. leeft 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 notice- ably outperform those employing solely Euclidean and exponential-map layers (see App. G.1.3 for an extensive analysis on SPD layers and App. G.1.5 for training times). Fig. E-middlle 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 rreviously-observe performance trend (see also App. G.I.I). Moreover, the geometric LNN trained over 8 integration steps outperforms those trained accelerations, and results in the most accurate predictions (see also App- G.1.4). A our results validate the effectiveness of considering additional inductive bias given by (7) the intrin- sic 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. G.1.2
Overall,
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 M (q: 07) = (gExp o (q) with Lr,sy_ 0.

---

Published as a conference paper at ICLR 2025
+ Cholesky
Cholesky By Cy
Expr Espr- Hum
ExPra Expr.
Expry GyroSpd...
Expp, CyroSpd,,
Expp, GymAl+RaFFig
Husi
10a
qll/ll
10-
10
10
10-1
100
10-1
11
100
7
10
0.0
0.5
1.0 time
1.5
2.0
0.5
1.0 time [s]
1.5
2.0
Draiil
Figure 2: 2-DoF pendulum: Left: Median acceleration prediction error for different LNNs and training set sizes Punin over 10 test trajectories. Shaded regions represent first and third quartiles. position and velocity predictions with respect prediction horizon Hust
Middle, right: Relative error of aumerically-interraa to the ground-truth trajectory
2000.
posiiio jaa
velccity (tad/s)
Dx2
D DoF
0.0
0.
0.2
0.3
0.4
0.7 0.1 0.2
0.5 0.5
[g]
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 updated with new initial condition (4u- do) every 0.025 (Husa 25).
5.2.1 COUPLED PENDULUM (16 DoFs)
As first ROM example, we consider 16-DoF pendulum. We ensure that the dynamics are ducible by constraining the motion the last 12-DoF as nonlinear combinations of the first 4-DoF. As such, the digensiimnsinaal pendulum is reducible to 4 latent dimensions.
Learning High-Dimensiona. 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 and the exponen- tial map basepoint to I (see App. Fo2 [ for implementation details). We consider models trained on acceleration data via (12). and with multi-step integration using (14) with Hmin = 8. Fig. 3 shows the AE reconstruction and RO-LNN predictions (Hes 25) for selected DoF's of Average testing errors for different Hust are reported in App. G.2. 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. Morcover, we compare our RO-LNN with full-order LNNs that directly learn high-dimensional dynamic parameters. As shown in Table even the best-performing FOM produces prediction errors that are orders of mag- nitude higher than the RO-LNN errors. Such erroneous acceleration predictions did not allow us to obtain stable velocity and position predictions. We also compare our models with the Lagrangian operator inference (L-OpInf) (Sharma & Kramer (024), and variant of L-OpInf that uses the same linear projection method with latent LNN for parameter identification (Sharma al. 2024). For fair comparisons, we replace the latent network used in (Sharma et al.] 2024) by a geometric LNN trained with multi-step integration (see App. F.2.). L-OpInf shows significantly higher errors than the RO-LNNs. While L-OpInf with LNN reaches lower errors, they still remain higher than for the similarly-traie trained RO-LNN (ODE), thus showcasing the limitations of linear subspaces.
trajectory.
Joint training. it 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
Table 1: Comparison of mean and standard deviation of prediction enrors over 10 test trajectories.
1.Oplff
1-OpInf (LNN, ODE)
RO-LNN (acc)
RO-INN (ODE)
(1.97 1.49) 10
(5.87 = 1.71) 10
(3.53 2.17) 10-1 13244=3.43) 10
(..99) 10
19, <l/110
(1.31 + 0.85) x 10 (2.45 + 0.96) x 10
(3.07 = 1.96) 10-1
0.39 x 10
8

---

Published as a conference paper at ICLR 2025
| gi1
9
9
v |INN,
FINNN
0.05
1.50
0:30
11.5
12.5
12.5
0.04
0:25
1.00
0:20
11.0 0,-
0.75
0.15
0.2
0:50
0:10
10.5
0:25
0.00
0.00
Figure 4: Median and quartiles of the errors of jointly-trained (on- Non-biorthogonal-manee manifoll vs overparametrized ()) and sequentially-trainee (I) models. Note that the latter have their own y-- axis in the left three plots, which is at least order of magnitude higher.
t 0.000s
= 0.750s
1.500s
t 2.250s
3.000s
= 3. 3.750s
-0.1
03
Figure 5: Predicted rope position (-) and ground truth (-) at selected timesteps for prediction horizon Hes 25. The grey circle depicts the circular trajectory of the end of the rope.
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 (/) ROM which sequentially trains the independent constrained AE and latent LNN, and (2) jointly-trained ROM that employs the overparamettrizaiion of (Otto (et al. 2023), and thus overlooks the geometry of biorthogonal AE weights. All models are trained with d via multi-ster integration with Huain = 8 and At
10 10-3s.
Fig. shows AE reconstruction errors (position and velocity) along with the high-dimensional and latent prediction eπrоrs averaged over 10 testing trajectories. We observe that the jointly- trained models outperform the sequentially-trainet one. Although the latter learns accurate dynamics 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 all loss components compared the RO-LNN trained the biorthogonal manifold. The overparametrized AE is also noticeably harder to train jointly with the latent LNN: In our experiments, 4 out of 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. G.2.2 for further comparisons).
5.2.2 ROPE (192 DoFs)
Next, consider high-dimensional deformable systems and learn the dynamics of 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. F.3 for details on the data and implementation). Fig shows the predicted rope configurations for a horizon of 25). Our model accurately predicts the high-dimensional dynamics of the rope (see also the predictions for selected DoFs in App. G.3 Notice that our attempts to train a full- order LNN were hindered by the high dimensionality of the system. Similarly, our attempts train O-Op did not result in any stable solution, thus substantiating the limitations of linear methods.
0.025 (Hhes
Latent space dimension. Selecting an appropriate latent dimension an important consideration for MOR. Here, we study the influence of the latent dimension on the performance of our RO-LNN. Table 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 nighhdimensional dynamics. However, the errors increase beyond certain latent dimension, indi- cating 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.

---

Published as a conference paper at ICLR 2025
Table Rope reconstruction errors and prediction loss components from (4)) for RO-LNNs.
4
1 q/
l
All
FNNNA
2.95 10 +324 10-2
4.63 10
10
3.06-228.00 10-2
3:00 10 1.34 10-2
6
8.19
10 +.11
10-1
1.52 x 10 2.92
10-2
2.40 10
4.57
10-2
1.81 10 5.69 x 10
10
4.75
x
10 5.733
10
6.53 x 10 1.62
10
1.43 10 3.45
10
2.87 10
8.92
14
6.16 10
6 10
5.68 x 10 8.48
10
9.22 10
1.19 10
5.06 10 1.08
10-1
s.00s
0.05
0.10s
0.15s
0.20s
0.25'
25 steps ground truth
full hol-honzzo
Figure 6: Predicted cloth configuration at selected times for 25-steps and 2500-steps horizons.
5.2.3 CLOTH (600 DoFs)
As a final example, we learn the dynamics of an even nigher-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 = 10 via multi-step integration on a set of 'Duainl = 8000 samples (see App. F.4 for details data and implementation). Fig. Gmiddle depicts the predicted cloth configurations for a horizon of s (Hees 25), showing that our model accurately predicts the cloth high-dimensional dynamics (see also the predictions for selected DoFs and errors in App. G.4 ppp
Full-horizon predictions. We evaluate the prediction quality of our model over a longer time hori- our model to predict the motion of the cloth for Hes 2500 timesteps. Fig. 6+boltom shows the cloth configurations predicted over the full horizon (see App. G.4 for predictions of selected DoFss and energy). Our model accurately predicts the cloth high- -igh-dinnnsoonal dynamics over the full horizon, suggesting that it successfully learned the constraints and actuation effects of the cloth falling on sphere.
zon. Namely, given an initial state (4) "
DISCUSSION
This paper proposed novel geometric network architecture for learning the dynamics of high- dimensional Lagrangian systems. By leveraging geometry and physies inductive bias, model pyysically-consistenn and conserves the reduced energy, infers interpretable reduced dynamic pa- rameters, 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 ge- ometry in our model. First, considering the geometry of Lagrangian systems allows nonlinear embeddings that ensure the preservation of the energy along solution trajectories and the interpretabiliit 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, addition to guaranteeing that they belong to the manifold of interest.
define
It is important to note that the prediction quality can only be 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 crucial, but non-trivial due to the unknown degree coupling in the FOM. In future work, we plan to extend our approach to applications beyond system identification. In particular, 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.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
This work was supported by ERC AdV grant BIRD, Knut and Alice Wallenberg Foundation and Swedish Research Council.
REFERENCES
Pitore-Antoine Absil, Robert Mahony, and Rodolphe Sepulchre.
on Matrix Manifolds. Princeton University Press, 2007. URL htttp:://press.pinnccttt edu/
absil
Joshua Barnett and Charbel Farhat. Quadratic approximation manifold for mitigating the Kol- mogorov barrier in nonlinear projection-bere based model order reduction. Journal of Computational Physics, 464, 2022. doi: https:///oi.org/ /0. 10.1001//..epp.... 1016//.jcp.202.. 11122.11144..
Gary Becigneul and Octavian-Eugen Ganea. Riemannian adaptive optimization methods. Intl. Conf. on Learning Representations (ICLR), 2019. URL htt.tttss//ppeeeeiieee net/
In
[ooumm Larum?dd==reiqq10K:
Rajendra Bhatia. Positive Definite Matrices. Princeton University Press, 2007. ISBN 0691 0691168533.
Aleksandar Botev, Andrew Jaegle, Peter Wirnsberger, Daniel Hennes, and Irina Higgins. Which
priors matter? Benchmarking models for learning latent dynamics. In Proceedings of the Neural
Information Processing Systems Track on Datasets and Benchmarks, volume 2021. URL
cappaper loteer__fiil files/
pdf
Nicolas Boumal. An introduction to optimization on smooth manifolds. Cambridge University Press,
2023. URL http:ttt://wwt /twt:::/www..r www.w:icw.was..m.amcommmmmmm www.ii.eet//ooooo
Patrick Buchfink, Silke Glas, and Bernard Haasdonk. Symplectic model reduction of Hamiltonian systems nonlinear manifolds and approximation with weakly symplectic autoencoder. SIAM
Journal of Scientific Computing, 45: 45:A289-A311, 2023. doi: 10. 1 11666657.
Patrick Buchfink, Silke Glas, Bernard Haasdonk, and Benjamin Unger. Model reduction on man-
ifolds: A differential geometric framework. arXiv preprint arXiv:2312.01963, 2024. URL
https://arxiv.cr wwwp:iii.rrgaaaaooo 22912019663
Kevin Carlberg, Ray Tuminaro, and Paul Boggs. Preserving Lagrangian structure in nonlinear model reduction with application structural dynamies. SIAM Journal on Scientific Computing, 37(2):
B153-B188 88, 2015. doi: 10. 00.1137/14095900
Elena Celledoni, Andrea Leone, Davide Murari, and Brynjulf Owren. Learning Hamiltonians of constrained mechanical systems. Journal of Computational and Applied Mathematics, 417: 114608, 2023. doi: 10.1016/j.cam. 2022. 114608.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural dinary differential equations. Neural Information Processing Systems (NeurIPS), vol- 31, 2018. URL httpsss/pproceeengg..nott neurips..........o ctcpaapp.ccppppppe wttpaper_flles/ppaeeeeeeeeeeeeee Files/paper/
Paper.
Miles D. Cranmer, Sam Greydanus, Stephan Hoyer, Peter W. Battaglia, David N. Spergel, and Shirley Ho. Lagrangian neural networks. ICLR Deep Differential Equations Workshop, 2020.
URL htR.htts:://axxtt
www.xiii.orrooooo.oo...... abs/2003. 04630
Thai Duong and Nikolay Atanasov. Hamiltonian-based neural ODE networks on the SE(3) manifold for dynamics learning and control. Robotics: Science and Systems (R:SS), 2021. doi: 10.15607/
ESS.202I.XVII.O
Giulio Evangelisti and Sandra Hirche. Physically consistent learning of conservative Lagrangian systems with Gaussin processes. In IEEÉ Conference on Decision and Control (CDC), pp. 4078- 4085, 2022. doi: 10.1 1099CCDDC50099222229991111

---

Published as a conference paper at ICLR 2025
Charbel Farhat, Todd Chapman, and Philip Avery. Structure preserrvng,, stability, and accuracy properties of the energy-conen conggy-consrrving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models. Intl. Journal for Numerical Methods in Engineering, 102, 2015. doi: (0.1002/mme.482)
Aasa Feragen and Andrea Fuster. Geometries and interpolations for symmetric positive defi- nite matrices, 85-113. Mathematics and Visualization. Springer, 2017. doi: 10.1007/ 978-3-319-61358-11
Mare Finzi, Ke Alexander Wang, and Andrew G. Wilson nian and Lagrangian neural networks via explicit constraints. In Neural Informa- fion Processing Systems (NeurIPS), volume 33, pp. 13880-13889. 2020. URL
Simplifying Hamilto
ortppp/ceee.aoeee.eeoro....o.oaa pdf
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural net- works. In Neural Information Processing Systems (NeurlPS), volume 32, 2019. URL
cwtppppppppppppppppppppppppppppe http://///pppp////////////////// Elles//paper2009/ 2bcoecadcooerddccccoo77oooo---pe
Jan S. Hesthaven, Cecilia Pagliantini, and Nicolò Ripamonti. Strecture-peeserving model order reduction of Hamiltonian systems. Intl. Congress of Mathematicians (ICM), volume 7, pp-
5072-5097, 2022. doi: 10.4171//CM202//111
Andreas Hochlehnert, Alexander Terenin, Steindór Samundsson, and Mare Deisenroth. Learning contact dynamics using physically structured neural networks. In Intl. Conf. Artificial Intel- ligence and Statistic (AISTATS), pp. 2152-2160, 2021. URL http::/ppocceediiicc//
mlr.
press//u3//hohhlleeetttttt/
Zhiwu Huang and Luc Van Gool. A Riemannian network for SPD matrix learning. In AAAI
Conf. on Artificial Intelligence, 2017. URL/tttps://oss.oea ierpaiiorrgiinndd hhppppppppaaa/a AAAI/
atticle/view/100655//1
Max Kochurov, Rasul Karimov, and Serge Kozlukov. Geoopt: Riemannian optimization in PyTorch. arXiv:2005.02879, 2020.
Kookjin Lee and Kevin Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. Journal of Computational Physics, 404:108973, 2020. doi: 10. jcpp://jjpp.21111111 1089939333..
Marco Lepri, Davide Bacciu, and Cosimo Della Santina. Neural autoencoder-basse structure- preserving model order reduction and control design for high-dimensional physical systems. IEEE Control Systems Letters, 8: 133-138, 2024. doi: 1..1000.SSSS.22233344222
Zhenhua Lin. Riemannian geometry of symmetric positive definite matrices via Cholesky decom- position. SIAM Journal on Matrix Analysis and Applications, 40(4):1353-1370, 2019. doi: 10.1137/18M1221000
Federico López, Beatrice Pozzetti, Steve Trettel, Michael Strube, and Anna Wienhard. Vector-valued distance and gyrocalculus the space of symmetric positive definite matri-
ces. In Neural Information Processing Systems (NeurIPS), volume 34, pp. 18350-18366. 2021. URL https::/proceedte wwurpps...........oo...........m cc wwtppppr.ppppp.eeeee www.s/ppper//e/////
F
Michael Lutter and Jan Peters. Combining physics and deep learning to learn continuous-time dynamics models. Intl. Journal of Robotics Research, 42(3):83-107, 2023. doi: 10.1177/
33783649231169499
Michael Lutter, Christian Ritter, and Jan Peters. Deep Lagrangian networks: Using physics as model prior for deep learning. Intl. Conf. on Learning Representations (ICLR), 2019. URL https:://penreviier hwt.:/..e/pppl... net/pdf Hp cc
12

---

Published as a conference paper at ICLR 2025
Xuan Son Nguyen. The gyro-structure of some matrix manifolds. In Neural Informa- tion Processing Systems (NeurIPS), volume 35, pp.
Curran Associates, Inc.,
2022. URL https://proeeed.gg www.:eeewdddnn....nne.e.e.....ee cc cc/paper Files/paper/2222. fi
Papeer ooneaaareeeetttttaooooooooeeeoar
Xuan Son Nguyen, Shuo Yang, and Aymeric Histace. Matrix manifold neural networks++.
In The Twelfth International Conference on Learning Representations, 2024. URL https:
7/ w.w.neeeie..eet............e http://ww.....oonmooo2----33333.
Samuel E. Otto, Gregory R. Macchio, and Clarence W. Rowley. Learning nonlinear projections
for reduced-order modeling of dynamical systems using constrained autoencoders. Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(11), 2023. doi:
Xavier Pennec, Pierre Fillard, and Nicholas Ayache. A Riemannian framework for tensor computing. International Journal of Computer Vision, 1):4416666, 2006. doi: 10070/ 1263-005-3222--.
Wilhelmus H. A. Schilders, Henk A. Van der Vorst, and Joost Rommes. Model Order Reduction: Theory, Research Aspects and Applications, volume of Mathematics in Industry. Springer Verlag, 2008. doi: 10.
Harsh Sharma and Boris Kramer. Preserving lagrangian structure in data-driven reduced-order mod- cling scrge-scale dynamical systems. Physica D: Nonlinear Phenomena, 462: 13412,, 2024. doi: https://doi.rr//00 wttpg/0011116000 1t..1//...yysd.2224444...
physd.2024.134125
Harsh Sharma, Hongliang Mu, Patrick Buchfink, Rudy Geelen, Silke Glas, and Boris Kramer. Sym- plectic model reduction of Hamiltomian systems using data-driven quadratic manifolds. Computer Methods in Applied Mechanics and Engineering, 417:116402, 2023, doi: 10...26... 116402.
Harsh Sharma, David A. Najera-Flores, Michael D. Todd, and Boris Kramer. Lagrangian operator inference enhanced with structure-preserving machine learning for nonintrusive model reduction of mechanical systems. Computer Methods in Applied Mechanics and Engineering, 423:116865 2024. doi: https://doi.ogg/1..001166. 10tp://...011.....2222444,, 16865.
Ryohei Shimizu, YUSUKE Mukuta, and Tatsuya Harada. Hyperbolic neural networks++. In Intl. Conf. Learning Representations (ICLR), 2021. URL hwLLttttes//eeeneviee...ee
Forum?id Ec85b0t
Yusuke Tanaka, Tomoharu Iwata, and Naonori Ueda. Symplectic spectrum Gaussin processes: Learning Hamiltonians from noisy and sparse data.
Neural Infor-
mation Processing Systems (NeurIPS), volume 35, pp. 20795-20808, 2022.
URL
https://ppoeedanss..eerss......p
25200406662266000622199922666663 Paper wwwfereno.eeeepceeeeeeeee
Maxime Thieffry, Alexandre Kruszewski, Christian Duriez, and Thierry-Marie Guerra. Control design for soft robots based on reduced-order model. IEEE Robotics and Automation Letters,
(I), 2019. doi: 10.1109//R..2018 28766744.
Emanuel Todorov, Tom Erez, and Yuval Tassa. Mujoco: A physics engine for model-based control. In JEFF/RSJ Intl. Conf. Intelligent Robots and Systems (IROS), pp. 5026-5033, 2012. doi: 200.///0...002..6366666.6
Yaofeng Desmond Zhong and Naomi Leonard. Unsupervised learning of Lagrangian dynamics from images for prediction and control. In Neural Information Processing Systems (NeurIPS), vol-
ume 33, pp. 10741-10752, 2020. URL httpps::/proceeediiini neur neur ps. www.ips.ccppppp...e.ee fil fttp://pprr/2000///////3ee3eeees
Yaofeng Desmond Zhong. Biswadip Dey, and Amit Chakraborty. Extending Lagrangian and Hamiltonian neural networks with differentiable contact models. In Neural In-
formation Processing Systems (NeurIPS), volume 34, pp. 21910-21922, 2021. URL
httpss///proceednngg..e..rr wwwring....uuuro.. httpp//errrilleeeppp//////////// b78888645913308833555111///ooof8 Pacsff--pperr.....
13

---

Published as a conference paper at ICLR 2025
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 S = {Σ E Sym" Σ > O} with Sym" denotes the set of symmetric R=xa matrices. The set S forms a smooth manifold of di- mension dim(S",) m(n+1), which can be represented as the interior of a convex cone embedded 2 in Sym". The tangent space TES+, at any point Σ E S is identified with the symmetric space Sym" R with origin at Σ.
The SPD manifold can be endowed with various Riemannian metrics. In this paper, we consider the affine-invariant metric (Pennec et al. (2006), which is widely used due to its excellent theoretical properties. The affine-invariant metric defines the inner product g TESt+ TyS", R given, for two matrices Tú, T2 E TzST +++ .....as as
( T2)x =
1).
(15)
The affine-invrriant metric defines the following geodesic distance
(A) log Σ AE -
3)ll.-
(16)
To map back and forth between the Euclidean tangent space and the manifold, we employ the cor- responding exponential and logarithmic maps Exps T=S",, - ST and Logz : ST, TySh
The exponential and logarithmic maps are then computed closed form as
A = Expx (S) = Σd exp Σ SE E
(17)
S = Logg (A) Σ log(2 Σ 3 トミ トコン
(18)
where exp(-) and log(-) denote the matrix exponential and logarithm functions. The parallel trans- port PTE-A TySt is defined as
TAS allowing us move elements between different tangent spaces
i T= PT (T) A T E-
(19)
with AE-A = - The exponential map. logarithmic map, and parallel transport are key operations for optimizing parameters on the SPD manifold in Sections 3 and
A.2 THE BIORTHOGONAL MANIFOLD
Pairs of matrices matrices &. V E Baxd with n 2 d > 1 that fullfil V T + I form the smooth manifold B,,d = {( v E Ruxd x Rex T = Ia] Otto et al. 2023 The biorthogonal man-
ifold B, is an embedded submanifold of the Euclidean product space Roxd "xd of dimension dim (B,,a) 2nd d² The tangent space at point (Φ), V)) E B, is given by
T v Brddd a {(V. W) E Baxd Raxd W T 0)
(20)
Vector can projected from the embedding space onto the tangent space via the projection opera- tion p Runxd x Raad
+ T( v Bu.d is given by
Proj
X Y) (x WAY 2
(21)
where is the solution of the Sylvester equation A A + TV) A - T
Riemannian operations the biorthogonal manifold are generally difficult obtain. Therefore,
leverage a retraction R x) T( 2) Ba,d - Budd as first-order approximation for the exponential
map. The retraction is formulated as
R x) (V, W) = -( (++ [ T + = v + W
(22)
A first-order approximation of the parallel transport operation the biorthogonal manifold PT [T (1) : T, Br.d - B is then defined via the successive applica-
tion of the retraction and projection as PT (.) Proj v2) R()...))
14

---

Published as a conference paper at ICLR 2025
B DRAWBACKS OF THE CHOLESKY DECOMPOSITION IN LNNS
The Cholesky decomposition M = LL with lower-triangular matrix L E R is commonly- employed to parametrize symmetric positive-definie matrices M S For instance, DeLaN Lut- ter et al.] (2019): Lutter & Peters (2023) enforces the symmetric positive-definitenss of the mass- inertia matrix M by inferring the lower-triangular matrix via a Cholesky network L(q; 0T) and reconstructing the mass-inertia matrix as M LL" While inferring L guarantees the symmetric positive-defimitenes of M, this method defines the distances between two SPD matrices based on the Euclidean distance between their Cholesky decompositions, therefore neglecting the curvature of
S This choice of distance suffers from the problematic swelling effect (Feragen & Fuster 2017]
Lin
2019), where the volume of the SPD matrices grows significantly while interpolating between
two
SPD matrices of identical volumes. Solutions inferred via a Cholesky decomposition will suffer
from this swelling effect and lead inaccurate predictions of the dynamies. These undesired effect can be avoided by directly inferring the mass-inertia matrix M in the SPD manifold ST equipped with the affine-invariant metric (Pennec et al.) 2006). To do so, we parametrize the mass-inertia matrix with a SPD network M (q; 0r) as explained in Section 3.1 and by optimizing its parameters using Riemannian optimization methods, as explained in Section 3.2
C LAYER TYPES IN SPD-NETWORKS
Here, we briefly introduce the different layers used in the SPD networks of Section
Euclidean Layers. Our SPD network leverages classical fully-comnected layers model functions that return elements the tangent space of manifold. The output of the I-th Euclidean layer x()) is given by
x (1) = (AA + bi),
(23)
with Ai E R and by E R the weight matrix and bias of the layer I, and nonlinear activation function of choice.
Exponential Map Layers. In the case of SPD manifolds, the exponential map layer is used to map layer inputs (1-1) E Sym" from the tangent space to the manifold. The layer output is given by
x Exp
(24)
with P E S", denoting the basepoint of the considered tangent space. One option is define P as equal to the identity matrix 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 local approximation to the manifold's curvature, the basepoint can instead be learned as network parameter P E 0r. s7 to potentially increase expressivity.
Gyrocalculus-basee FC SPD Layers (GyroAl). Together with gyrocalculus formulations for the SPD manifold, López et al. (202T) and Nguyen (2022) presented a version of fully connected layers for SPD data. With @ as addition and @ as scaling operations for SPD gyrovectors.
x = (Ai & x + B
(25)
defines the 1-th gyrolayer output (1) St_ for the input (1-1)) E S Here, the gyrovector (X)) with * denoting pointwise multiplication and Matrices A,, B: ST.
operations correspond to A @ X = exp( A
X B = B
part of the network parameters ST.
Gyrospace hyperplane-based FC SPD Layers (GyroSpd..). Inspired by their equivalent in hy- perbolic spaces (Shimizu et al. 2021), another version of a gyrocalculus- FC layer for SPD networks was proposed by Nguyen al. (2024). 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 from specific point V, that is obtained from the layer input X (1-1) S via the gyrocalculus formulation of the hyperplane equation
hyperplane that is located in the origin of the manifold
v(X
(eBIe) (1-1)) AI),
(26)
15

---

Published as a conference paper at ICLR 2025
that is parametrized by the network parameters A Bi S
The additive inverse fulfills
ex = x- The final layer output is given by
exp( (xp(((((-1)) * I,); (
(27)
where (x())) = /2 for + j and Iy(s)) for j.
ReEig Nonlinearities. The ReEig layers were first introduced by Huang & Gool (2017), and were extended to more general context in Nguyen (2022). Even though 8 already classifies linear space, this type of layer adds the option of adding further nonlinearities to the network. The layer output is defined
1 = spu X (1)-1)) -ddia(maxx(( V))UT,
(88)
with eigenvalue decomposition of (1-I) Udiag (V)ia.(V)U The constant acts rectification threshold to avoid smaller eigenvalues and hereby regularize the output. Under the inveriinvariant distance, this nonlinearity layer is analogous to ReLu nonlinearity in the Euclidean network layers.
D CONSTRAINED AE: IMPLEMENTATION DETAILS
The constrained AE is implemented following OOtto al., 2023). To guarantee the projection prop- erties, the nonlinear activation functions employed in the encoder and decoder network and a+ must satisfy a_ o 0+ = id. To do so, they are defined by
bzz V2
1
V2 sin(a) Sin(aa)cosaa)) cos(a)
2x
(zi)
-
+ 2a,
(29)
asin(a)
with
fa csc2 sec? (a),
(30)
b csc + sec?
The activations then resemble smooth, rotation rotation-symmetric versions the common leaky ReLu ac- tivations. The parameter 0 a < sets the slope of the activation functions. Throughout experiments, we set a =
For each layer the weight matrices given by the pair (11, 1 Bu (see App. A.2 In the proposed RO-LNN, we train these weights via Riemannain optimization on the biorthogonal manifold, analogous the Riemannian optimization the SPD manifold presented in Section 3.2 Notice that the optimization of the AE parameters = = { V bi takes place on product of Euclidean and biorthogonal manifolds.
When learning dynamic latent space for Lagrangian MOR, the AE should Iearn the lifted embedding (8 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. second derivative of the encoder and decoder networks enables training on acceleration data. Full-dimensional acceleration data is encoded via the reduction map Po as q = dp q + d &ppoliq and latent accelerations are a = de qu + d' (q.4) q. This requires the expensive
reconstructed via the embedding fe q computation of second-order derivatives, but is nevertheless part of our experiments.
Given dataset of T observations {qi, qi, q. T.)TT the AE is trained on loss
LAAE
q.. q.112 llq. 4.11 0. 11.
(31)
of position reconstructions q. velocity reconstructions q. and acceleration reconstruction ğ. 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. 2023), where the biorthogonality of weights is achieved
16

---

Published as a conference paper at ICLR 2025
92 depla d docka
q
P2 dpolg
T(q. q
4p
q
q
F(
q
V(()
depe4 4
Figure 7: Flowchart the forward dynamics of the proposed reduced-order LNN trained on ac- celeration via (12). The reduction mappings and embeddings of the Lagrangian ROM are depicted in blue and parametrized via constrained AE with biorthogonal layers. The ROM dynamics are learned via latent geometric LNN depicted in orange. The mass-inertia matrix is parametrized via the SPD network described in Section 3
via the introduction additional penalty losses. In this case, the overall loss only reaches a min- imum value if the overparametrized weight matrices actually fulfill the biorthogonality conditions for each layer I. The overparametrized constrained AE loss is given by
LAAE up LoAE LOAE, my
with LAAE. ry
+ (1)
$ O
I
||
(32)
with the loss LeAE, reg regularizing the biorthogonality of the AE weights. Note that this over- parametrization does not guarantee the biorthogonality condition, in contrast to the Riemannian approach. In practice, the additional regularization loss feAE, reg required an additional learning rate, resulting in more possible combinations of hyperparameters.
TRAINING ON ACCELERATION
Fig. illustrates the complete forward dynamics of the proposed RO-LNN trained acceleration by minimizing the loss 12
F ADDITIONAL EXPERIMENTAL DETAILS
This section presents additional details on the experimental setup of Section Position and velocity predictions obtained via Euler forward integration.
F.1
2-DoF PENDULUM OF SECTION 5.1
F.1.1 SIMULATION AND DATA GENERATION
For this experiment, we use 2-DoF pendulum implemented in MUJOCO (Todorov et al. 2012) The pendulum connected via hinge joints. Both links are modeled as cylinders of radius r; 0.025 m, length 4 = 0.4 m, and mass m, = 0.1 kg. The simulated environment is entirely free of dissipative effects. We generate two different datasets D Dpnnd = {fq qo.k: 4n.ks Tn, k== consisting of N trajectories with K samples per trajectory. The state evolution is simulated in MUJOCO with an RK4-solver with timestep Al = 10-3, s over a time interval I I= [0, 7]s.
1K IN
randomly sampled from the in- terval g: = 0) E [0, 303 for each DoF = {1, 2). The initial velocities are set to q(t = 0) = 0. Each trajectory recorded for T 2s. All models presented in Section 5. in the main text trained on this dataset.
Dataset 1: Unactuated pendulum. The initial configurations
Sine-Tracking pendulum. Each joint is controlled to follow reference trajectory via an inverse dynamics torque control law. The reference trajectories are set to follow sinusoidal shapes as Gref, A, s2n fit + Φ), with Φ; aresin 20.0 A.
Dataset 2:
4:,0 E [0,30] o A, E [1,30]* %
fo E [+, s11, and i = (1,2}. Cier.s and qher.: are obtained as first and second derivatives of Weet with respect time. Initial positions and velocity
match the reference trajectories as
17

---

Published as a conference paper at ICLR 2025
q.(t = 0) st 0) and (t 0) = quet, (t 0) Each trajectory is recorded for T 3.5s.
Results of models trained on this datasets are presented in App. G.1
F.1.2 ARCHITECTURES TRAINING
Geometric LNN architecture. For all models, the potential energy network (qq: 0y) consists of Lv 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 M (q: 0r), which differ in their components JR, Nisps and gs.+ as follows:
OR: The Euclidean component consists of Lr, Euclidean layers on the tangent space of the SPD manifold with SoftPlus activations and 64 neurons per layer.
Iixpi The basepoint S", is either as the origin learned as a network parameter Po E S
952,: The SPD component consists of LT,s._ hidden SPD layers. We consider two dif- ferent FC SPD) layers, namely GyroAI (25).
2.. and GyroSpd,+ (26). Both FC layers can be
augmented with ReEig nonlinearity (88). The regularization constant of the nonlinear-
ity is uniformly set as = 1 x 10-4. Note that for the SPD layers, the common notion of number of neurons per layer layer width not directly applicable. As the layers as transformation S"; + S via network parameters 0r,sT they can be interpreted carrying single transformation per layer, or as transformations of manifold dimension dim(S m(++1)
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: (2)
The kinetic and potential energy networks are independent in (Lutter & Peters 2023). Both
networks consist of 2 Euclidean SoftPlus layers of 64 neurons. Hereby, the potential energy network fully equivalent to the ones employed in the geometric LNN. The kinetic energy network returns a Cholesky decomposition L(q) of the mass-inertia matrix M (2) The kinetic and potential energy networks share parameters, i.e., Or ( Oy as in (Lutter al. 2019). The MLP consists of 2 hidden Fuclidean SoftPlus layers of 64 neurons, and separate output layers return the potential energy and the Cholesky decomposition.
Training on Acceleration Data. Each training dataset Duain consist of randomly-sampled data from 30 trajectories from one of the datasets presented in Section F.I.I The models trained by minimizing the LNN loss (6) for 3000 epochs with architecture-spenne learning rates. If not further specified, both gR and gST, consist of each 2 layers.
Training on multi-step Integration. For these experiments, we randomly sample recordings form 30 trajectories, which then form the training dataset D {19)10; to.- {q)) to. [[]) Note that this training dataset depends on the ODE solver's timestep At and desired prediction horizon H. The models trained by minimizing the loss
jAt 11âp:(tes)) ( (((t)) (1g.)(ti.s) gi (1.,)))) + 1101 (33)
LNN HN
for 1500 epochs with architecture-specffic learning rates. Notice that the loss (33) is the equivalent of (14) 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 I.
Long-term prediction Experiments. In those experiments, use both models trained celeration data and with multi-step integration. In both scenario, the training datasets contain 'Dmiin 8000 samples and are trained for 2000 epochs.
18

---

Published as a conference paper at ICLR 2025
Table Nonlinear combinations for DoFs 16 of the coupled pendulum.
DoF
f(a:: 92. 43 Q4)
46
23 cos(
91 + 0.1 (12(
97
q4 cos
48
+ 43 (q2
99
Q10
99190
f11
sin(
Q12
0.44334
f13
-0.9q1
92 + qs
2q
Q14
-3 -3sin
Q15
-293
Q16
-0.9
F.2 COUPLED PENDULUM (16 DoFS) OF SECTION 5.2.1
F.2.1 SIMULATION AND DATA GENERATION
In this experiment, we consider 16-DoF coupled pendulum. The pendulum is connected via hinge joints. All links are modeled as capsules of radius r; = 0.05 m, cylinder length 0.5 m, and mass my 1 kg. The simulated environment is entirely free of dissipative effects. We generate a dataset D Ppndi6 = {1qn,k Kn.k. 4r.k. Ta.k ] === 1n=1 consisting of N trajectories with K samples per trajectory. The first 4-DoF are simulated in MCJOCO (Todorov et al. 2012) with an RK4-solver with timestep At = 10 10-3, s over a time interval == [0. 3|s. while the last 12 DoFs are constrained nonlinear combinations of the first as described in Table The initial configurations of the first 4 DoFs are randomly sampled from the interval qt = 0) E [0,30] o. The initial velocities are d(t = 0) 0.
K L.V
F.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 = n0 pendulum dynamics. The constrained AE is implemented with 4 biorthogonal encoder and decoder layers Po
match the known underlying dimension of the high-dimensional
(2) R R and o I R of sizes nı [8, 16, 16, 16). For the latent geometric LNN, both the potential and kinetic energy network consist of 2 hidden Fuclidean layers of 64 neurons. The kinetic energy network employs exponential map layer with basepoint Expr-
We consider two versions of the RO-LNN. The first model is trained on acceleration data by mini- mizing loss 12). The second model trained with multi-step integration by minimizing (14). We consider Humin 8 latent integration steps and an integration time constant of 10-3 s. The losses are optimized using Riemannian Adam. For the first model, we use learning rate of x 10 for the AE parameters E. 1 x 10 for the LNN parameters 0, and regularization = 1 x 10 6 For the second model, the learning rate for the parameters of the LNN Both models are trained until convergence, i.e., for 4000 and 3000 epochs, respectively.
x 10 with 7 2 x 10-5
We compare our RO-LNN with (/) a full-order geometric LNNs, (2) L-OOpnff (Sharma Kramer 2024). and (3) L-Opinf with latent LNN, as described in (Sharma et al.] 2024). For the full- order geometric LNNs, we only report results of the best-performing FOM out kinetic energy network consists of LT.R Euclidean hidden SoftPlus layers with 128 neurons, an exponential map layer with basepoint P I, and no SPD layers, i.c., Lr,s, = The potential energy network is likewise composed of 2 Euclidean SoftPlus layers 128 neurons. The full-order LNN trained for 3000 epochs. The learning rate for the LNN parameters is set to 1 x 10 and = x 10
models. The
consider projection onto 4-dimensional linear sub- space, as for our RO-LNN. We implemented L-OpInf in Python and solved the optimization problem leading to the reduced equations of motion with CVXPY. As described in (Sharma & Kramer 2024
For L-OpInf (harma & Kramer 2024).
19

---

Published as a conference paper at ICLR 2025
set the low-dimensional mass-inertia matrix as constant identity matrix, i.e., M Id. In our experiments, we optimize for 15000 samples of ground-truth position data.
For L-OpInf with latent LNN (Sharma al., 2024), we considered the same projection dimensional linear subspace as for the L-Oplni bascline. For fair comparisons, we enhance the latent LNN compared to that used in (Sharma al. 2024). In our experiment, the latent dynamic parameters are learned by geometric LNN trained via the multi-step integration loss (14), for which
4-
use the same hyperparameters for the RO-LNN trained with the same loss. Note that, for both versions of L-OpInf, the model is fed with ground truth data every Hes 25 steps to maintain consistency comparison with the RO-LNN. We used the same Euler forward integration scheme as with our other models.
When investigating joint training Section 5.2 the previously-desenine RO-LNN trained with multi-step integration is compared to (/) version trained via the overparametrization of biorthogo- nal weights, and (2) a model that was trained separately, by sequentially training the AE and latent LNN. The overparametrized model uses the overparametrized loss (32) with fAAE, reg multplied by weighting constant x 10 For the sequentially- -ereiinnell-ttrineedd model, the first 3000 epochs were trained using only the AE components fA the loss (14). For the consecutive 3000 epochs, only the LNN components LnN were active of (14). All models were trained on |aiainl 24000 samples, which were randomly sampled from 20 trajectories. The testing dataset consists of 10 trajectories.
F.3 ROPE (192 DOFS) OF SECTION 5.2.2
F.3.1 SIMULATION AND DATA GENERATION
The rope is implemented in MUJOCO via elastic cable. Over total length of m, the rope consists of i = 64 equally-spaced capsule-shaped masses m, 0.1 kg and r; 0.02 m. Twisting and bending stiffness are set to 5 x 10S 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 all directions of the workspace. For training and testing, we generate the datasets Drope {(qm,k. Qn,k: trajectories with K samples per trajectory. The generalized coordinates of the system are chosen to be the Cartesian positions qi (are, Vi, 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 This mimicks a handheld manipulation task.
1.W consisting of N
For each reference trajectory, a radius ry is sampled randomly from r; E [0.05, 0.4] m. To achieve non-planar motions, the circles are tilted by random angle 0 E [-90, 90] The end of the rope is set to track such a reference at an angular velocity 120 as via PD control. The state evolution simulated with an RK4-solver with timestep At = 10 - s over a time interval I = [0, 3.8|s.
F.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 latent ge- ometric LNN. The constrained AE is implemented with 4 biorthogonal encoder and decoder layers (1) R + RT and 49 Ru=-1 + Rn of sizes n4 [64, 64, 128, 192]. We consider latent P0
(1)
spaces of dimension d [4, 6, 10, 14). In all scenarios, the kinetic energy network of the latent geometric LNN consists of hidden Fuclidean layers of 64 neurons with SoftPlus activations, and exponential map layer Exp, The potential energy network also consists of 2 hidden Euclidean layers of 64 neurons with SoftPlus activations.
All models are trained on Hnnn 8 integration steps by minimizing the loss [4 with Riemannian Adam (Becigneul Ganea 2019). All models are trained with a learning rate of x 10 for the AE parameters E, 1 x 10 for the LNN parameters 0, with regularization 7 1 x 10 1 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, record 10 full trajectories.
20

---

Published as a conference paper at ICLR 2025
F.4 CLOTH (600 DoFs) OF SECTION 5.2.3
F.4.1 SIMULATION AND DATA GENERATION
The deformable thin cloth is modelled in MUJOCO flexible composite object with 200 masses of m, = 0.1 kg equally spaced over a width of 0.1 m and length of 0.2 m. For this model, the generalized coordinates are given by the Cartesian positions 4 s. Mis 2,] of each masses' center of mass in the world frame.
For training and testing, generate datasets Peloth Km kn T= ] consisting of N = 20 trajectories during training, and N = 10 trajectories during testing. The state evolution simulated with an RK4-solver with timestep At 10-4 s over a time interval resulting in K' 2500 samples. The trajectories contain recordings of the cloth falling on sphere from a height of 0.13 in the center above the origin of the sphere. To vary scenarios, the radius of the sphere is andooll-samplee from E [0.02, 0.12) m.
[0, 0. 0.25|s,
F.4.2 ARCHITECTURES & TRAINING
The RO-LNN is composed of a constrained AE with layer sizes n4 [32, 64, 128, 600] for the = 1,.., 4 encoder and decoder layers Po
(1) R + R= and o R - R and latent space dimension no 10. The kinetic energy network of the latent geometrie LNN consists of 2 hidden Euclidean layers of 64 neurons with SoftPlus activations, and an exponential map layer Expy. The potential energy network also consists of 2 hidden Euclidean layers of 64 neurons with SoftPlus activations.
The featured model is trained on Hrain integration steps via the loss 14 optimized using Riemannian Adam with learning rate of 5 10 for the AE parameters E, and 2 x 10 for the LNN parameters 0, and regularization
x 10 THe model is trained for 3000 epochs.
G ADDITIONAL EXPERIMENTAL RESULTS
G.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
G.1.1 ADDITIONAL RESULTS ON THE EXPERIMENT OF SECTION
Trajectory prediction. Fig. complements Fig. midddee right by depicting the long-term position and velocity predictions of the different geometric LNNSs considered 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.
Energy conservation. Fig. displays the total energy of the different LNNs of Figs. 2 along the predicted trajectory, As the pendulum is unactuated, the energy levels - that can be learned up to constant should be constant over the predicted trajectory. We observe that all LNNs approx- imately conserve the total predicted energy E, with the DeLaNs showing the highest fluctuations. Note that, since the predicted trajectories based on Euler forward integration, we do not expect the total energy of the system to be perfectly preserved.
G.1.2 SINE-TRACKING 2-DoF PENDULUM
We reproduce the experiment of Section 5.1 2-DoF pendulum dataset introduced in App. F1..
different dataset, namely the sine-tracking
Fig. 10 depicts the acceleration prediction errors for selected architectures. Similarly for the un- actuated dataset, the LNNs with Cholesky networks are outperformed by the geometric LNNs, and the low-data regime. In contrast to the unactuated pendulum dataset, the geometric
21

---

Published as a conference paper at ICLR 2025
Cholesky
Cholesky Or n Bv Expp, GyroAI+ReEig position [rad]
Expy Expr.
Expro Exps Hinnn
Expn, GyroSpd
Expr, , +ReEig
Hum
Ground truth
velocity rad/s]
-0.5
0.05
0.0
0.00
-0.5
1.75
tine
lime [s]
Figure 8: Numerically-integraee position and velocity predictions and ground truth trajectory for LNNs trained on the unactuated 2-DoF pendulum dataset.
Cholesky
Cholesky
Expx Expr
ExP ExP1: Huin
Exppy GyroSpd
Expp, GyroSpd.. +ReBig
Expa, GyroAI+ReBig
Huzi
Ground truth
5.20
9.5
0.55
0.00
-0.0
-000
Line
[+]
Figure 9: Total energy & (up to a constant) of different LNNs.
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 obtain accurate trajectory prediction for the sine-tracking dataset. Generally, the results on this second dataset are aligned with those pre- sented in Section 5.1 and validate the effectiveness of considering the geometry of the mass-inertin matrix and training the models via multi-step integration when learning Lagrangian dynamics.
G.1.3 INFLUENCE OF SPD LAYERS
Next, we evaluate the influence of SPD layers within the kinetic energy network M (q: 0r) of the geometric LNN. We consider the network architectures described in App. and train the models with acceleration data. We vary the type of SPD layers, as well the number of hidden layers LT,RR Ir,sz, 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 accel- eration prediction errors for the unactuated and sine-tracking datasets are given in Tables and respectively.
Despite the increased performance of geometric LNNs over DeLaNs, the reported prediction errors do not allow to observe any clear influence of the choice of SPD layers, exponential map base- points, 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 GyroAl, and ReEig layers in SPD network may be due practical issues related to the training of these
22

---

Published as a conference paper at ICLR 2025
Cholesky
Cholesky Or n Bv
Expy Expr.
Expro Exp2. Hunn
Expr, GyroSpd...
Expr, GyroSpd,, EReEig
Expr, GyroAI+ReEig
Hum
109
1000
2000
3000
4000
5000
6000
7000
[irriil
Figure 10: Median acceleration prediction error for different LNNs and training set sizes Dimin over 10 testing trajectories. The models were trained on the sine-tracking 2-DoF pendulum dataset. Shaded regions represent first and third quartiles.
Cholesky
Cholesky O2 li By Expp, GymAI+ReFig
Expz Expr
Exp Expr: Husi
Exppy- GyroSpd...
Expp, GymSpd,, +Refig
Ruis
Ground truth
al/ll
100
10
10
10
2
3
2
3
time
[
time
Figure 11: Relative emrors of numerically-integrabee position and velocity predictions w.r.t. ground truth trajectory over prediction horizon Hest for LNNs trained on the sine-tracking 2-DoF pendulum dataset.
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, indicated by 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 replace unstable autodiff computations, endowing the SPD manifold with a different Riemannian metric.
Due to the limited performance improvements and increased computational complexity (see App. G.1.5) of the SPD layers), we consider SPD networks M (q: 0T) (gEsp o (R) (q) composed of Euclidean and exponential map components when learning reduced-order Lagrangian dynamics in the experiments of Section 5.2
Table Mean and standard deviation of the acceleration prediction errors 11 El for geometric LNNs with different SPD network architectures trained on various training sizes on the unactuated 2-DoF pendulum dataset.
DDnial
1000
2000
3000
5000
7000
Exp,
LI 1. Lr.80 3
5.95 5.69
5,22 3.78
1.83 2.61
1.86 2.98
0.92 = 1.31
GyroSpd,,,
LT LT.
19.5
3.66
4.43 2.96
2.71 3.22
1.52 1: 1.69
L1A 1,1,,88
1.18
1:17
1.75
Expry
L1A Ln. =2
L.37 2.66
3.3423772
2.08 3.02
0.89 1.11
1.44 2.02
IΛ 1, =3
4.00 + 2.0
1.73
1.77 2.400
1.00 = 2.58
1.10 1.42
GyroSpd.,
LT ITSS. =2
1.66 2,81
1.72 2.19
1.80 2.52
1.78 2.21
0.84e + 1.25
Expre
Ix 11,
2.85
L In.5
1.99 1.994.229 2.23
2.38
2.62
1.77
1.14 1.31
ExPr
2,31 1 2.33
7,28 1. 8.43
1.52
2.08
1.32 1.91
1.13
GyroAI+RoEig
L1 LT =2
1.77 + 2.34
Ban
1,99
=: 2.70
nan
1.30
2.17
Expry
L7x LT,80
2.76
12.49
L.1
21.2
19.3
1.28 1.00
1.07 1.45
LTx LT,
23

---

Published as a conference paper at ICLR 2025
Table 5: Mean and standard deviation of the acceleration prediction errors و6 - اll for different geometric LNNs with different SPD network architectures trained on various training sizes on the sine-tracking 2-DoF pendulum dataset.
Puuia 3000
1000
2000
5000
7000
Expr
Lia = 11,8, =3
3.73 + 1.46
122 = 1.12
3.95 2 1.113
3.37
3.51
2.05 E 2.41
GyroSpd...
Lr LT =2
10.9 8.62
3.95
4.50
3,89
4.49
2.80 2.68
2.18 1: 2.29
Expra
TTa 1,1,, 1T.8 las,
12.0
3.58 4.58
3.95
3.44
4.84
3.66 + 4.40
2.75 2.80
1.49 4 1.47
ExPr
1,1, m3
4.59 1.00:44000 4.58
4.17 4.43
4.15 4.448
71.0
79.6
6.80 5.23
, +RoEig
LA = las. =2
13.5 15.3
1.05
4.26
3.95 + 4.38
5.78
4.91
1.27 + 1.25
Expr
Lra = 1,,1,
1:29 1.33
1.55 = 1.17
4.13::1:11
1.71 1.88
3.31 3.015
LT LT,8,
2.82
L1 L1.5% Lra Less, =2
6.73
1.00 1.45
3.33 4.30
6.74
Expr
9.71
3.79
1.11
3.67 + 1.60
2.36 2.58
19.1 58.2
Espra
L1a = 11. =3
1.02 - 4.06
4.11
2.35
2,08
1,86 1:113
LT = 1Ts5 =2
3,99 01 4.66
3.40
4.55
2.78
2.68
1.68 1. 1.75
Table 6: Mean and standard deviation of the Hist 25-steps velocity prediction errors for a geo- metric LNN trained via multi-step integration with different horizons Hiain and datasets.
Dataset
Hean =1
Hean
Huam
Huin 12
Unactuated
4.61 x 10-1 = 6.83 10-1
6.64 x 10-2 1.73 x 10-1
4.80
10-2 = 1.17 x 10"1
4.80 x 10-2 1.13 10-1
Sinc-tracking
10 == 1.07
10
8.87 x 10 2.05
1.25
10
1.89
10
1.36
10
2.22 10-2
1.1.4 INFLUENCE OF MULTI-STEP INTEGRATION DURING TRAINING
We assess the influence the number Husann of steps considered when training our geometric LNNs. The featured networks trained as described in F.1.2 For this experiment, the architectures
as described for multi-step integration in section F.1.2 TJrpoorres the velocity prediction errors averaged over 10 testing trajectories for models trained different prediction horizons Hmin 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 Huain. Ultimately, we select Hurain by trading off the increase in performance and the training time. The later is investigated next.
G.1.5 RUNTIMES
Table shows the average training times of several LNNs trained on datasets of |Duain = 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_,++R +ReFig layers as they are the computationally most demanding SPD layer of these layers.
combination. The featured model consists of sT
Importantly, the well-performing geometric LNN with identity basepoint and LT.s = 0 trains almost fast as DeLaN with its Cholesky layers. Instead, integrating CyroSpd,.+ReEinn 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 well-performing LNN with Euclidean layers and the Exp, 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
Exp,
Expp, 121
Expp, GyroSpd ReHig
Expy, Hyain 8
92
189
105
24

---

Published as a conference paper at ICLR 2025
Table Mean and standard deviations the reconstructed and latent prediction errors for different prediction horizons with RO-LNNs trained on acceleration and via multi-step integration.
Σ (n(t;) ((t))
fNNN
Σ ANNN
Has
ODE
acc
1.86 x 10 + x 10-4
2.65 x 10 + 3.18 x 10
4.06 x 10 8 x 10-2
1.96 x 11" + #2144 10"
1.70 x 10 + 4.94 x 10
9.22 + 3.76 x 10-3
Mun
25
ODE
acc
2.05 x 10 4.277 10-4
x =13.12 x
1.66 x 10 100-14..80 4.85 x 10-1
2.25 x 10"50 10"
x 4 10
1.45 - 4.28
Hes
50
acc
+3100
10
10"
10
2
ODE
4.03 10 8.62
10
5.41
10 1.70 109
2.76 x x
1111100--111110 1 10
G.2 ADDITIONAL, RESULTS ON THE 16-DOF COUPLED PENDULUM EXPERIMENT OF SECTION 5.2.1
This section presents additional results on learning the high-dimensional dynamics of 16-DoF pendulum with geometric RO-LNN.
G.2.1 ACCELERATION VS MULTI-STEP TRAINING
Here, we further evaluate the impact the losses used to train geometric RO-LNNs. Table g shows the average acceleration errors of trajectory predictions of RO-LNNs trained on accelerations via multi-step integration with horizon Haain = 8. We observe that the RO-LNN trained via multi-step integration consistently outperforms the RO-LNN trained on accelerations for all prediction horizons Heat 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 traiming 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.
G.2.2 CONSTRAINED STANDARD AUTOENCODERS
We further investigate the effectiveness our geometric constrained AE whose parameters are obtained via Riemannian optimization compared to two baselines, namely (I) constrained AE op- timized via the overparametrizatioo of (Otto et al. 2023), and (2) a regular (unconstrained) AE. All AEs are trained to reduce set of position data Duain {q from 15 trajectories of the 16-DoF pendulum to d 4-dimensional latent space. The constrained AEs are implemented as in Section 1 and described in App. F2.. The regular AE consists of 4 Fuclidean layers with ReLu activations of sizes nn = [64, 64, 64, 16) for 1 1, 1,.,," 4 encoder and decoder layers (1)
Pg R + R and 40 R + R Our constrained AE is trained using Riemannian Adam with learning rate 5 10-2, The overparametrized and regular AEs are trained using Adam with learning rate 5 x 10 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 con- strained AEs consistently outperform the regular one. Morcover, the Riemannian optimization leads to lower reconstruction errors than the overparametrivaation, especially in low-data regimes.
G.2.3 COMPUTATIONAL SPEEDUP OF ROM W.R.T. FOM
Despite that we consider scenarios where the FOM is unknown and that trainings of a full-order LNN were unsuccessful (see Se2tion(5.2.I1, we aim at providing idea of the computational effort of our model compared to the evaluation of the FOMs. To do so, we symbolically derive the La- grangian equations of motion with known physical quantities of the considered 16-DoF pendulum from Section F.2. to obtain FOM. We compare the wall-clock time of evaluation of this FOM against the ODE-trained RO-LNN, For both models, we consider a roll-out of 3000 timesteps from the same initial conditions qf 0s) and q(( 0s) via Euler forward integration with timestep At 30 Averaged over 10 runs on the same local CPU, we achieve evaluation times of 113.59 for the FOM, and 1.57 for the ROM. We conclude that the ROMs obtained with our pro- posed method not only enable structure- structure-preserrvee preserving learning of high-dimensionall system dynamics, but also significantly improve computational efficiency compared to evaluating FOMs.
25

---

Published as a conference paper at ICLR 2025
0.5
0.3
0.0
1000
2000
3000
4000
5000
6000
7000
train dataset size
Figure 12: Median and quartiles of position reconstruction error for the regular (-), the biorthogonal constrained AE (-), and the overparametrized constrained AE (-) for RO-LNN trained on different sizes training sets.
position [m]
velocity [m/s]
8.02.
0.00
Dop
-20000
DoP 02
0.02
0,05
0.00
DoP 33
-0.29
0.15
0.0
DoP 121
0.00
-0.05
-0.05
Dof 155
0.035
0.0
0.000
DoF 192
0.05
01.00
00.25
0.50
0.25
1.00
1.75
0.00
0.75
1.00
1.25
1.50
175
time [g]
time [8]
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 (qo- qo) every 0.025 (Hess = 25).
G.3 ADDITIONAL. RESULTS ON THE 192-DoF ROPE EXPERIMENT OF SECTION 5.2.2
Fig. 13 shows the AE reconstruction and the RO-LNN predictions (Hest = 25) for selected DoFs of test trajectory. We observe that the the RO-LNN accurately models the nigh-dimensional dynamics of the rope.
G.4 ADDITIONAL RESULTS ON THE 600-DoF CLOTH EXPERIMENT OF SECTION 5.2.3
Table shows the reconstruction and prediction errors of our model averaged over 10 testing trajec- tories. Moreover, Figs. 14 and show the AE reconstruction and the predictions of the RO-LNN of selected DoFs of test trajector for horizons Hist 25 and Hist 3000, respectively. Notice that the ground truth only contains actuation data for the first 0..55 (i.e., 2500 timesteps). Therefore,
26

---

Published as a conference paper at ICLR 2025
Table Cloth reconstruction and Hust 25-steps prediction errors over 10 testing trajectories for the RO-LNN with d = 10 trained via multi-step integration with the loss 14
19 1--9
ol
2.47 x 10
1,50 10" 1.81
10" 1.51 10" 1.81 x 10" 5.07 10
1.00 10
posiiio [rad]
velccity (rad/s]
0.0025
DoF
0.0000
-0.1
0.0025
DoF 200
0.0025
0.0000
0.0
000025
Dof 300
0.0
-0.1
2000
0.0085
0.05
Dop
0.0000
0.00
7-0.0025
0.00
Dof 5A0
0.0025
0.0000
--0.0025
DoF 800
0.0
-0.1
0.05
0.20
time [
time
Figure 14: Cloth position and velocity predictions (-) from the RO-LNNs trained via multi-step integration and ground truth (-). The corresponding AE reconstructions (....) depicted for completeness. The model is updated with a new initial condition (go. go) every 0.025 (Hest = 25).
reasonably assume that the cloth remains still for the additional time horizon laug (0.2 00.25, 0.30 and set the ground truth states and torques equal to their last recorded values. Fig. 16 disp 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 com- plex high-dimensional deformable systems.
27

---

Published as a conference paper at ICLR 2025
posation [ra]]
velocity [rad/s]
0.005
ToP
-0.005
0.000
-0.1
0.005
DoP
0,000
0.0
-0.005
0.0
DoF.
-001
0.000
DoF
10,000
0.005
DoF 500
0.005
0.1
0,000
DoPP n60
0.0
-001
Line
Sinze
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 (qu- qo) for 0.3 HHs 3000). The grey-shaded areas indicate the interval during which ground truth data extrapolated from the last observation.
0.05
5. 0.00
-0.05
0.00
0.05
0.10
0.15 time [s]
0.20
0.25
Figure 16: Predicted (-) and ground truth (-) latent energy E (up to a constant) over time. The grey-shaded areas indicate the interval during which ground-trun ttrth data are extrapolated from the last observation.
28


Reviews:
[{"review_id": "iXJmjVeLSd", "reviewer": "Reviewer_4zWX", "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.", "questions": "Most of this reviewer's questions is dealt with in the weakness section above.\n\nAdditionally, how could we connect the satisfaction of the Euler-Lagrange equations of the latent dynamics to that of the original dynamics?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "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.", "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.\n\nThis identification of the contribution should be explained clearly, and the performance of RO-LNN should be compared to the prior works in the experiments.\n\nThere 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.\n\nFor 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.\n\nWhile 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.\n\nFor 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.\n\nIn 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.\n\nFor 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "kQQjGodpHZ", "reviewer": "Reviewer_7VkD", "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.", "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?\n2. 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?\n3. (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?\n4. 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?\n5. 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?\n6. How was the time step $\\Delta t=10^{-3}$ chosen? How would choosing larger time step affect the roll-out prediction error?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The paper presents experimental results on increasing complexity of the case\n2. It shows better results compared to Lagrangian Neural Network", "weaknesses": "Major:\n1. 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.\n2. 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.\n3. 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.\n\nMinor comments:\n1. Please avoid using the unscientific notation like $1e^{-5}$. Instead, please use notations like $1\\times 10^{-5}$.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "44SlvIoduw", "reviewer": "Reviewer_rYLi", "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.", "questions": "see the weakness part.", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 1, "strengths": "The strengths are listed as followed\n- Utilizes geometry and physics as inductive biases, leading to a physically consistent model to improve the interpretability. \n- Effectively learns reduced-order dynamics for high-dimensional rigid and deformable systems with accurate long-term predictions.\n- Joint optimization of latent space and reduced dynamic parameters enhances performance, particularly in low-data regimes.\n- Incorporates intrinsic geometry of the problem, aiding in the preservation of energy and interpretability of the model.\n- Outperforms standard autoencoders by employing a constrained auto-encoder trained with Riemannian optimization.", "weaknesses": "### Contribution and Originality\n\n- 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. \n\n- 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].\n\n- 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.\n\n\n### Comparison\n\n- 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.\n\n### Computational Complexity\n\n- 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. \n\n\n[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.\n\n[2] Buchfink, Patrick, et al. \"Model reduction on manifolds: A differential geometric framework.\" Physica D: Nonlinear Phenomena 468 (2024): 134299.\n\n[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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "BSlEbPDIod", "reviewer": "Reviewer_wYEc", "summary": "The paper proposes a reduced-order learning model with the preservation of the state geometry for high-dimensional Lagrangian systems.", "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.\n\nQ2. The background and design motivation should be emphasized for readers unfamiliar with Lagrangian systems. See my Weakness points 1, 2.\n\nQ3. The writing logic should be improved. See my Weakness points 5, 6.\n\nQ4. The computational cost should be illustrated and reported. See my Weakness point 8. \n\nQ5. The result presentation can be improved. See my Weakness points 9 and 10.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "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.", "weaknesses": "1.\tIn 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.\n\n2.\tIs 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? \n\n3.\tWhat’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}$?\n\n4.\tWhat is $V(q;\\theta_V)$ in Line 162? Since there are massive notations, the author is advised to complete the proofreading. \n\n5.\tIt’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?\n\n6.\tThe 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. \n\n7.\tNotations in Eq. (8) are unclear to the reader. What are $l$, $h$, $x$, and PT? Can the author align them with the previous notations? \n\n8.\tIt 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.  \n\n9.\tThe legends in Fig. 2 are hard to align with the curves. Please refine this figure. \n\n10.\tWhy is the error in Table 1 very high? I am not sure if this is a validated result.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)