Neural Contractive Dynamical Systems
Authors: Hadi Beik Mohammadi
Søren Hauberg
Georgios Arvanitidis
Nadia Figueroa
Gerhard Neumann
Leonel Rozo
Venue: ICLR 2024 spotlight (2024)
Keywords: learning from demonstration
dynamical systems
contraction theory
Primary area: applications to robotics, autonomy, planning
TL;DR: The paper introduces the Neural Contractive Dynamical System (NCDS), providing a flexible and stable approach for learning contractive dynamics, and extends it to high-dimensional systems with obstacle avoidance capabilities.

Abstract: Stability guarantees are crucial when ensuring that a fully autonomous robot does not take undesirable or potentially harmful actions. Unfortunately, global stability guarantees are hard to provide in dynamical systems learned from data, especially when the learned dynamics are governed by neural networks. We propose a novel methodology to learn \emph{neural contractive dynamical systems}, where our neural architecture ensures contraction, and hence, global stability. To efficiently scale the method to high-dimensional dynamical systems, we develop a variant of the variational autoencoder that learns dynamics in a low-dimensional latent representation space while retaining contractive stability after decoding. We further extend our approach to learning contractive systems on the Lie group of rotations to account for full-pose end-effector dynamic motions. The result is the first highly flexible learning architecture that provides contractive stability guarantees with capability to perform obstacle avoidance. Empirically, we demonstrate that our approach encodes the desired dynamics more accurately than the current state-of-the-art, which provides less strong stability guarantees.

Published as a conference paper at ICLR 2024
NEURAL CONTRACTIVE DYNAMICAL SYSTEMS
Hadi Beik-Mohammadi Søren Hauberg', Georgios Arvanitidis', Nadia Figueroa',
Gerhard Neumann', and Leonel Rozo!
1 Bosch Center for Artificial Intelligence (BCAI), 2 Karlsruhe Institute of Technology (KIT),
3 Technical University of Denmark (DTU), Uniiersity of Pennsylvania (UPenn).
Emails:
madii.biik-mol meikmmoohmmaaii, Iconel. dee.de.bssee bosch.coc, [sohau, gearl @duu.dkk nadiafig @ seas, upunn.euu, gerhard. .eeuaannneemanneee @ kit.edu
ABSTRACT
Stability guarantees are crucial when ensuring that fully autonomous robot does not take undesirable or potentially harmful actions. Unfortunately, global stability guarantees hard to provide in dynamical systems learned from data, especially when the learned dynamics are governed by neural networks. We propose a novel methodology to learn neural contractive dynamical systems, where neural architecture ensures contraction, and hence, global stability. To efficiently scale the method nigh-dimensional dynamical systems, develop a variant of the variational autoencoder that learns dynamics in a low-dimensional latent representation space while retaining contractive stability after decoding. We further extend our approach to learning contractive systems the Lie group of rotations to account for full-pose end-effector dynamic motions. The result the first highly flexible learning architecture that provides contractive stability guarantees with capability to perform obstacle avoidance. Empirically, we demonstrate that our approach encodes the desired dynamics more accurately than the current state-of- the-art, which provides less strong stability guarantees.
STABILITY IN ROBOT LEARNING
Deploying fully autonomous robot requires guarantees of stability to ensure that the robot does not perform un- wanted, potentially dangerous, actions. Consider the robot in Fig. 1, which is trained from demonstrations (black dots) to execute a sinusoidal motion (orange). If the robot is pushed away (red perturbation) from its planned tra- jectory, then should still be guaranteed to reproduce a motion similar to the demonstration pattern (blue). Man- ually designed robot movements could achieve such sta- bility guarantees, but even skilled engineers struggle hand-code highly dynamic motions (Billard et al., 2016). In contrast, learning robot dynamics from demonstrations has shown to be an efficient and intuitive approach for encoding highly dynamic motions into robot's repertoire (Schaal et al., 2003). Unfortunately, these learning based approaches often struggle to ensure stability as they rely on the machine learning model to extrapolate in a con- trolled manner. In particular, it has turned out to surprisingly difficult to control the extrapolating behavior of neural network models (Xu et al., 2021), which hampers stability guarantees.
baaseng-bbased
Figure Robot motion executed via a neural contractive dynamical system (NCDS).
Stability is commonly ensured through asymptotic or contraction guarantees. Asymptotic stability ensures that all motions converge to fixed point (known as the system's is attractor) (Khansari-Zadeh & Billard, 2011; 2014; Neumann & Steil, 2015; Zhang et al., 2022b;a). This is suitable when the only requirement is that the robot eventually reaches certain configuration, e.g. its end-effector at a specific goal position. Many tasks, however, require the robot to dynamically follow desired trajectories, e.g. in flexible manufacturing. human-robot interaction, in entertainment settings.

---

Published as a conference paper at ICLR 2024
In these cases, asymptotic stability is insufficient. A stronger notion of stability is provided by contraction theory (Lohmiller& & Slotine, 1998), which ensures that all the path integrals regardless of their initial state incrementally converge over time. Unfortunately. the mathematical requirements of a contractive system are difficult to ensure in popular neural network architectures.
Existing contractive learning methods focus low-capacity models fitted under contraction con- straints. Ravichandar et al. (2017) and Ravichandar & Dani (2019) both use Gaussian mixture models to encode the demonstrations and model the dynamical system with Gaussian mixture regression (GMR), which is fitted under contraction constraints. Blocher et al. (2017) also use GMR to model the system dynamics but instead include a stabilizing control input that guarantees a local contractive behavior. Likewise, Sindhwani et al. (2018) model the dynamical system with a well-crafted kernel, under which optimization remains convex w.r.t. contraction constraints. Note that the above works impose the contraction constraints during optimization, which complicates model training.
Developing high-capacity (deep) models difficult as off-the-shelf neural network architectures provide no stability guarantees. Instead, current approaches split the task into two steps: first, learn a non-contractive dynamical system, and then estimate Riemannian metric (Sec. 1). parametrized by a neural network, under which the system becomes contractive (Tsukamoto et al., 2021a; Dawson et al., 2023). The approach, thus, requires training two separate networks, which is impractical and comes with additional challenges. For example, Sun et al. (2020) learns the metric by regularizing towards a contractive system, such that stability can only be ensured near training data. Alternatively, Tsukamoto & Chung (2021) develop convex optimization problem for learning the metric from sampled optimal metries. Unfortunately, accessing such optimal metrics is a difficult problem in itself. Additionally, Kozachkov et al. (2022) introduce network of small and simple recurrent neural networks with built-in contractive behavior as an alternative to a single, larger network.
In this paper. we propose the neural contractive dynamical system (NCDS) which is the first neural network architecture that is guaranteed to be contractive for all parameter values. This implies that NCDS can easily be incorporated into existing pipelines to provide stability guarantees and then trained using standard unconstrained optimization. As demonstration data may be high dimensional, we secondly propose an injective variant of the variational autoencoder that allows learning low- dimensional latent dynamics, while ensuring contraction guarantees for the decoded dynamical system. Thirdly, as most complex robotic tasks involve full-pose end-effector movements (including orientation dynamics), extend approach to cover the Lie group 83(3). Finally, show that NCDS can trivially be modified to avoid (dynamic) obstacles without sacrificing stability guarantees.
2 LEARNING CONTRACTIVE VECTOR FIELDS
Our ambition is a flexible neural architecture, which is guaranteed to always be contractive. We, however, first introduce contraction theory this is prerequisite for our design.
2.1 BACKGROUND: CONTRACTIVE DYNAMICAL SYSTEMS
To begin with, assume autonomous dynamical system o; - f(x)). where Ø, RD is the state variable, f RR= + RDD is a C² function and é; = dæ/ar denotes temporal differentiation. As illustrated in Fig. 1, contraction stability ensures that all solution trajectories of a nonlinear system f incrementally converge regardless of initial conditions Dg. Zg. and temporary perturbations (Lohmiller & Slotine, 1998). The system stability can, thus, be analyzed differentially, i.c. two nearby trajectories converge to one another. Specifically, contraction theory defines a measure of distance between neighboring trajectories, known the contraction metric, in which the distance decreases exponentially over time (Jouffroy & I. Fossen, 2010: Tsukamoto et al., 2021b).
ask
Formally, an autonomous dynamical system yields the differential relation dæ (2)6z, where J J(x) = af/2æ is the system Jacobian and 8æ is a virtual displacement (i.e., an infinitesimal spatial displacement between the nearby trajectories a fixed time). Note that we have dropped the time index t to limit notational clutter. The rate of change of the corresponding infinitesimal squared distance dx 8x is then,
d
(8x 6xx) 28x dे dt
28x J())(x
(I)
follows that if the symmetric part of the Jacobian ((x) is negative definite, then the infinitesimal squared distance between neighboring trajectories decreases over time. This is formalized as:
2

---

Published as a conference paper at ICLR 2024
Definition 1 (Contraction stability (Lohmiller & Slotine, 1998)) An autonomous dynamical sys- tem Æ = f(z) exhibits a contractive behavior if its Jacobian J(x) = 61 of/0a aw is negative definite, or equivalently if its symmetric part is negative definite. This means that there exists constant 0 such that Ex dæ converges to zero exponentially at rate 2r. This can be summarized as,
BT S.t. Ya, [/2 + - O.
(2)
The above analysis can be generalized to account for more general notion of distance of the form Sæ M (z) where M(x) (x) is a positive-definite matrix known as the (Riemannian) contraction -
metric (Tsukamoto al., 2021b; Dawson et al., 2023). In our work, we learn a dynamical system f so that it is inherently contractive since its Jacobian [(x) fulfills the condition given in Eq. 2. Consequently, it is not necessary to specifically learn a contraction metric as corresponds to an identity matrix. For our purposes, we, thus, define contraction guarantees following Eq. 2.
2.2 NEURAL CONTRACTIVE DYNAMICAL SYSTEMS (NCDS)
From Definition 1, we seek a flexible neural network architecture, such that the symmetric part of its Jacobian is negative definite. We consider the dynamical system ø f(æ), where x R denotes the system's 's state and f IR + RD is parametrized by a neural network. how to impose a negative definiteness constraint on network's 's Jacobian without compromising its expressiveness. To achieve this, first design neural network J, representing the Jacobian of our final network. This produces matrix-valued negative definite outputs. The final neural network will then be formed by integrating the Jacobian network. Specifically, we define the Jacobian as,
not obvious
fz())= - -Jo (x) Jo(æ) + <IDp),
where Jo: IRD - Dx a neural network parameterized by 0, 6 TR+ is a small positive constant, and Tn is an identity matrix of size D. Intuitively, Jo can be interpreted as the (approximate) square root of J,. Clearly, Jy is negative definite as all eigenvalues are bounded from above by -6.
Next, we take inspiration from Lorraine & Hossain (2019) and integrate Jo to produce a function f. which is implicitly parametrized by 0, and has Jacobian Jy. The fundamental theorem of calculus for line integrals tells us that we can construct such function by line integral of the form
④ = f(x) = -000 '''' J, (c (t), č(x, wo)dt,
(4)
c(x, t, wo) = (1 - t) æ to and ö(t, xo. x) wo (5) where Dy and ato f (mr) represent the initial conditions of the state variable and its first-order time derivative, respectively. The input point æo can be chosen arbitrarily (e.g. the mean of the training data or can be learned), while the corresponding function value ite has to be estimated alongside the parameters 0. There- fore, given a set of demonstrations denoted as D = {x,, }. our objective is to learn set of parameters along with the initial conditions Ny and d.g, such that the integration in Eq. 4 enables accurate reconstruction of the velocities d; given the state x;,
The learned vector field (grey) and demonstrations (black). Yellow and green trajectories show path integrals starting from demonstration starting points and random points, respectively.
Figure 2:
The integral in Eq. 4 is similar to neural ordinary differential equations (Chen et al., 2018b), with the subtle difference that it is a second-order equation as the outcome pertains to the velocity state æ. We can, thus, view this system as a recond-order neural ordinary equation (Norcliffe et al., 2021) and solve it using off-the-shhe numerical integrators. The resulting function f will have a negative definite Jacobian for any choice of 0 and is consequently contractive by construction. In other words, we can control the extrapolation behavior of the neural network that parametrizes
our dynamical system f via the negative-deffnitenes of (x). We call this the neural contractive dynamical system (NCDS), and emphasize that this approach has two key benefits: (1) we can use any smooth neural network Jo unconstrained optimization, unlike previous approaches. Figure 2 shows an example NCDS vector field, which is clearly highly flexible even if it guarantees contractive stability.
base model, and (2) training then be realized using ordinary
3

---

Published as a conference paper at ICLR 2024
3 LEARNING LATENT CONTRACTIVE DYNAMICS
Learning highly nonlinear contractive dynamical systems in high-dimensional spaces is difficult. These systems may exhibit complex trajectories with intricate interdependencies among the system variables, making it challenging to capture the underlying dynamics. In the experiments, we show that NCDS works very well for low-dimensional problems, but when only limited data is available, the approach becomes brittle in higher dimensions. A common approach in such cases is to first reduce the data dimensionality, and work before in the resulting oi-diiensional latent space (Chen et al., 2018at 8a; Hung et al., 2022; Beik-Mooammadi et al., 2021). The main difficulty is that even if the latent dynamics are contractive, the associated high-dimensional dynamics need not be. This we solve next.
3.1 BACKGROUND: DEEP GENERATIVE MODELS
Variational autoencoders (VAEs). The main goal of deep generative models is to approximate the true underlying probability density p(æ) given a finite set of training data in an ambient space X, by considering a lower-dimensional latent space Z. In particular, the variational autoencoder (VAE) (Kingma & Welling, 2014; Rezende et al., 2014) is a latent variable model, often specified through a prior p(z) = N( 0, Ia) where Z, and likelihood P& (a) []) N(x [z)(2), I (=)) with X. Typically, the mean and the variance of the likelihood are parametrized using deep neural networks Hs 2 - X and o3 2 + E with parameters Ø, and Ip and Ia being identity matrices of size D and d, respectively. These neural networks are trained by maximizing the evidence lower bound (ELBO) (Kingma & Welling, 2014). The latent variable z is approximated using variational encoder (x) and is often interpreted as the low- now-dimensional representation of an observation æ. In our work, we use a VAE to provide low-dimensional representations of individual points along observed trajectories in order to learn a latent contractive dynamical system.
Injective flows. A limitation of VAES is that their marginal likelihood is intractable and we have to rely on a bound. When dX' = dim(Z), we can apply the sfange-of-variablee theorem to evaluate the marginal likelihood exactly, giving rise to normalizing flows (Tabak & Turner, 2013). This requires the decoder to be diffeomorphic, i.e. a smooth invertible function with smooth inverse. In order to extend this to the case where dim(X)= > dim(Z), Brehmer & Cranmer (2020) proposed an injective flow, which implements a zero-padding operation (see Sec.3.2 and Equation 6) on the latent variables alongside a diffeomorphic decoder, such that the resulting function is injective.
3.2 LATENT NEURAL CONTRACTIVE DYNAMICAL SYSTEMS
As briefly discussed above, we want to reduce data dimensionality with a VAE, but we further require that any latent contractive dynamical system remains contractive after it has been decoded into the data space. To do so, we can leverage the fact that contraction is invariant under coordinate changes (Manchester & Slotine, 2017; Kozachkov et al., 2023). This means that the transformation between the latent and data spaces may be generally achieved through a diffeomorphic mapping.
Theorem (Contraction invariance under diffeomorphisms (Manchester & Slotine, 2017)) Given a contractive dynamical system á f()) and a diffeomorphism vē applied on the state E D. the transformed system preserves contraction under the change of coordinates y = ((x). Ou Sz.
Equivalently, contraction is also guaranteed under a differential coordinate change by = Dw
Following Theorem 1, we learn a VAE with an injective decoder μ 2 - X. Letting M = (Z) denote the image of μ. then μ is diffeomorphism between 2 and M,, such that Theorem I applies. Geometrically, l spans a d-dimensional submanifold of X on which the dynamical system operates.
Here we leverage the zero-padding architecture from Brehmer & Cranmer (2020) for the decoder. Formally, an injective flow ; 2 + X learns an injective mapping between a low-dimensional latent space Z and higher-dimensionall data space X. Injectivity of the flow ensures that there are no singular points or self-intersections in the flow, which may compromise the stability of the system dynamics in the data space. The injective decoder μ is composed of zero- paro-padding operation on the latent variables followed by series of K invertible transformations 9k- This means that,
/ 9K o 91 Pad,
(6)
4

---

Published as a conference paper at ICLR 2024
where Pad(z) [21 ... 0...0]
0]' represents D-dimensional vector with additional D d' zeros. We emphasize that this decoder is an injective mapping between 2 and f((2 C X, such that a decoded contractive dynamical system remains contractive.
Specifically, we propose to learn a latent data representation using a VAE, where the decoder mean l4 follows the architecture in Eq. 6. Empirically, we have found that training stabilizes when 1 (x), [aaz( where Iie" ser1 is the
the variational encoder takes the form (=(=(x) N( me approximate inverse of lig given by,
Mg Unpad 91
9K
where Unpad:
1R + Rd removes the last D-d dimensions of its input as an approximation the inverse of the zero-padddee parddppdding operation. We emphasize that an exact inverse is not required to evaluate a lower bound of the model evidence.
It is important to note that the state a solely encodes the positional information of the system, disregarding the velocity æ. In order to decode the latent velocity & into the data space velocity á. we exploit the Jacobian matrix associated with the decoder mean function μgg, computed as Jug (z) = /8z/ This enables the decoding process formulated below,
a Jug (z)
(8)
The above tools let us learn a contractive dynamical system on the latent space Z, where the contraction guaranteed by employing the NCDS architecture (Sec. 2.2). Then, the latent velocities & given by such contractive dynamical system can be mapped to the data space X using Eq. 8. Assuming the initial robot configuration 26 is in M (i.e., In = (zo)), the subsequent motion follows contractive system along the manifold. If the initial configuration 2o is not in M, the encoder is used to approximate a projection onto M, producing 2y / (zx). Note that the movement from Zu to [=z) need not be contractive, but this is finite-time motion, after which the system is contractive.
3.3 LEARNING POSITION AND ORIENTATION DYNAMICS
So far, we have focused Euclidean robot states, but in practice, the end-effector motion also involves rotations, which do not have an Euclidean structure. We first review the group structure of rotation matrices and then extend NCDS to handle non- -on-Euclidean data using Theorem 1.
Orientation parameterization. Three-dimensional spatial orientations can be represented in several ways, including Euler angles, unit quaternions, and rotation matrices (Shuster, 1993). We focus on the latter. The set of rotation matrices forms Lie group, known as the special orthogonal group S0(3) {R R3x3 R R det R2) 1}. Every Lie group is associated with its Lie algebra, which represents the tangent space at its origin (Fig. 3). This Euclidean tangent space allows us operate with elements of the group via their projections on the Lie algebra (Solà et al., 2018). In the context of (3), its Lie algebra so(3) is the set of all 3 x 3 skew-symmetric matrices [r) x- This eeemersymettric matrix exhibits three degrees of freedom, which can be reparameterized as 3 3-dimensional vector
[r,, Ty, r2] E R3.
First coster.
We can map back and forth between the Lie group 8O(3) and its asso- ciated Lie algebra (3) using the logarithmic and exponential maps, denoted Log (3) + s0) and Exp- : s0(3) - 80(3). respectively. We provide explicit formulae for these in Appendix A.2. Due to wrapping (e.g. 360 rotation corresponds to 09), the exponential map is surjective. This implies that the inverse, i.e. Log, is multivalued, which complicates matters. However, for vectors r E so (3), both Exp and Log are diffeo- morphic if ||| (Falorsi et al., 2019; Urain al., 2022). This ball B, is known as the first cover of the Lie algebra and corresponds to the part where no wrapping occurs.
Exp
Log
SO(3)
Figure 3: Aspects of the Lic group 83(().
NCDS on Lie groups. Consider the situation where the system state a rotation, i.e. x E 80(3), which we seek to model with a latent NCDS. From a generative point of view, we first construct a decoder μ + s3(3) with outputs in the Lie algebra. We can then apply the exponential map to generate rotation matrix R E 80(3), such that the complete decoder becomes Exp on.
For training, the latent velocities simply estimated by numerical differentiation the latent
5

---

Published as a conference paper at ICLR 2024
Irlx
Encoder
Decoder
TC41
r++1x
Ju
Pr
Pri
80(3)
(3(3)
R,
R+t1
P
R
j
Re-1
Ja
U+++1
O
2g 10000
Figure 4: Architecture overview: single iteration of NCDS simultaneously generating position and orientation dynamics. (A) VAE (pink box): The encoder processes the concatenated position-orientation data pt, yielding resulting vector that subsequently divided into two components: the latent code z (yellow squares) and surplus (gray squares). The Unpad function in Equation 7 removes the unused segment. The unpadded latent code is fed to the contraction module and simultaneously padded with zeros (white squares) before being passed to the injective decoder. (B) Contraction (blue box): The Jacobian network output, given the latent codes, reshaped into square matrix and transformed into negative definite matrix using Equation 2. The numerical integral solver then computes the latent velocity I Later, using Fq. 8, á is mapped to the input-space velocity via the decoder's. s Jacobian Jug.
Unfortunately, even if is injective, we cannot ensure that Exp ou is also injective (since Exp is surjective), which then breaks the stability guarantees of NCDS, Here we leverage the result that Exp is diffeomorphism as long we choose a decoder architecture such that μ : Z. - B- is injective and have outputs on the first cover, then Exp ou is injective, and stability is ensured (see Appendix A.5). Note that during training, observed rotation matrices can be mapped directly into the first cover of so(3) using the logarithmic map. When the system state consists of both rotations and positions, simply decode to higher dimensional variables and apply exponential and logarithmic maps on the appropriate dimensions. Figure 4 depicts this general form of NCDS.
restrict ourselves to the first cover of so(3). Specifically, if
3.4 OBSTACLE AVOIDANCE VIA MATRIX MODULATION
In real-world scenarios, obstacle avoidance is critical to achieving safe autonomous robots. Thus, the learned dynamical system should effectively handle previously unseen obstacles, without interfering with the overall contracting behavior of the system. Fortunately, NCDS can easily be adapted to perform contraction- peretiion-reeevinngg.. obstacle avoidance using the dynamic modulation matrix G from Huber et al. (2022). This approach locally reshapes the learned vector field in the proximity of obstacles, while preserving contraction guarantees. To avoid obstacles effectively, it is imperative to know both the position and geometry of the obstacle. This makes obstacle avoidance on the VAE latent space difficult as we need to map both obstacle location and geometry to the latent space 2.
Instead, we directly apply the modulation matrix the data space X of the decoded dynamical system. Specifically, Huber et al. (2022) shows how to construct modulation matrix G' from the object location and geometry. such that the vector field & (( 2 is both contractive and steers around the obstacle. With this minor modification, NCDS can be adapted to avoid obstacles. provide the details in Appendix A.3. Note that this is tailored for obstacle avoidance the end- effector level. In order to extend this to multiple-limb obstacle avoidance, the modulation matrix must be refined to incorporate the robot body, and NCDS must trained using joint space trajectories.
4 EXPERIMENTAL RESULTS
To evaluate the efficiency of NCDS, we consider several synthetic and real tasks. Comparatively, we show that NCDS the only method to scale gracefully to higher dimensional problems, due to the latent structure. We further demonstrate the ability to build dynamic systems the Lie group of rotations and to avoid obstacles while ensuring stability. Neither of the baseline methods have such capabilities. Experimental details are in Appendix A.4, while ablation studies of NCDS modeling choices are in Appendix A.6.3. Furthermore, the videos are available at: http://siisssogll.eooooooe wwwp:/gwwww...o.eeeeeeeeeo www.m//wwww/eeaalcotttttttoooo/o
https://ites..oogge
Datasets. There currently established benchmarks for contraction-stabl robot motion learning, so we focus on two datasets. First, we test our approach on the LASA dataset (Lemme al., 2015), often used for benchmarking asymptotic stability. This consists of 26 different two-

---

Published as a conference paper at ICLR 2024
Figure Visualization of the LASA-2D dataset: Gray contours represent the learned vector field, black trajectories depict demonstrations, and orange/green trajectories illustrate path integrals starting from the initial points of the demonstrations and plot corners. The magenta circles indicate the initial points of the path integrals.
dimensional hand-written trajectories, which the robot is tasked follow. To evaluate method, we have selected 10 different trajectories from this dataset. To ensure consistency and comparability. preprocess all trajectories all stop at the same target state. Additionally, points of each demonstration, to ensure that the only state exhibiting zero velocity is the target state. To further complicate this easy task, task the robot follow 2, or stacked trajectories, resulting in 4 (LASA-4D), or 8-dimensional (LASA-8D) data, respectively. We further consider a dataset consisting of trajectories from 7-DoF Franka- Emika Panda robot using a Cartesian impedance controller.
omit the initial few
Metries. We use dynamic time warping distance (DTWD) as the established quantitative measure of reproduction accuracy w.r.. a demonstrated trajectory, assuming equal initial conditions (Ravichandar et al., 2017; Sindhwani et al., 2018). This is defined as,
DTWD(T Tax Tw')
Σ min (El(72)
(d(T Tat)
Σ min jEllrz)
JEI(T2))
(El(T.))
where Tx and Tu- are two trajectories (e.g. a path integral and a demonstration trajectory), d is distance function (e.g. Euclidean distance), and ((T) is the length of trajectory 7.
Baseline methods. Our work is the first contractive neural network architecture, so we cannot compare NCDS directly to methods with identical goals. Instead, we compare to existing methods that provide asymptotic stability guarantees. In particular, Euclideanizing flow (Rana et al., 2020), Imitation flow (Urain al., 2020) and SEDS (Khansari-Zadee) & Billard, 2011).
4.1 COMPARATIVE RESULTS
We first evaluate NCDS on two-dimensional trajectories for ease of visualization. Here data is sufficiently diw-dnmeniional that we do not consider latent structure. Figure shows the learned vector fields (gray contours) for 10 diverse trajectory shapes chosen according to their difficulty from the LASA dataset, covering a wide range of demonstration patterns (black curves) and dynamics. We observe that NCDS effectively captures and replicates the underlying dynamics. We compute addi- tional path integrals starting from outside the data support (green curves) to assess the generalization capability of NCDS beyond the observed demonstrations, Remarkably, NCDS successfully generates plausible trajectories even in regions not covered by the original demonstration data. Moreover, the yellow path integrals, starting from the initial points of the demonstrations, show that approach can also reproduce the demonstrated trajectories accurately. This is evidence that the contractive construction is viable approach to controlling the extrapolation properties of the neural network.
Table 1 shows average DTWD distances among five generated path integrals and five demonstration trajectories for both NCDS and baseline methods. For the two-dimensional problem, both Eu- clideanizing flow and Imitation flow outperform NCDS and SEDS, though all methods perform quite well. To investigate stability, Fig. displays the average distance over time among five path integrals starting from random nearby initial points for different methods. We observe that only for NCDS does this distance decrease monotonically, which indicates that it is the only contractive method.

---

Published as a conference paper at ICLR 2024
NCDS
EF
SEDS
J
m
Figure 7: Path integrals generated using different methods on LASA-8D. Black curves are training data, while yellow are the learned path integrals. To construct the 8D dataset, we have concatenated different 2D datasets.
To test how these approaches scale to higher-dimensional settings, consider the LASA-4D and LASA-8D datasets, where we train NCDS with two-dimensional latent representation. From Table 1, is evident that the baseline methods quickly deteriorate as the data dimension increases, and only NCDS gracefully scales to higher dimensional data. This is also evident from Fig. which shows the training data and reconstructions of difffernnt LASA-8D dimensions with different methods.
Mean distance
NCDS
IF
SEDS
200
400
600
800
Time
Figure 6: Average distance between random nearby trajectories time for LASA-2D. Only NCDS monotonically decreases, i.e. is the only contractive method.
These results clearly show the value of having a low-dimensional latent structure (more details in Appendix A.6.5).
To further compare our method, we consider a robotic setting where we evaluate all the methods on real dataset of joint-space trajectories collected on 7-DoF robot. The last column of Table 1 shows that only NCDS scales gracefully to high-dimensional joint-space data, and outperforms the baseline methods. The supplementary material includes video showcasing the simulation environment in which each robot executes a drawing task, depicting V-shape trajectories on the table surface.
4.2 OBSTACLE AVOIDANCE
Unlike baselines methods, our NCDS framework also provides dynamic collision-free motions. To assess the obstacle avoidance capabilities of our approach, we conducted a couple of experiments on the LASA dataset. The results of this experiment are reported in the third panel of Figure 8. As observed, the presence of an obstacle that completely blocks the demonstrations is represented by a red cirele. Remarkably, NCDS successfully reproduced safe trajectories by effectively avoiding the obstacle and ultimately reaching the target. Moreover, the last panel Figure shows the obstacle avoidance in action the robot's end-effector navigates around the orange cylinder. These results experimentally show how the dynamic modulation matrix approach leads to obstacle-free trajectories while still guaranteeing contraction.
4.3 ROBOT EXPERIMENT
To demonstrate the potential application of NCDS roboties, we conducted several experiments on a 7-DoF Franka-Emika robotic manipulator, where an operator kinesthetically teaches the robot drawing tasks. The learned dynamics are executed using Cartesian impedance controller. The details of the implementation and the specific network architecture are provided in Appendix A.4 and A.6.4.
LASA-2D
LASA-4D
LASA-8D
DoF robot
Euclideanizing flow
0.72 0.12
3.23 1 0.34
10.22 0.40
5.11 + 0.30
Imitation flow
0.80 + 0.24
0.79 + 0.22
4.69 + 0.52
2.63 + 0.37
SEDS
1.60 + 0.44
3.08 + 0.20
4.85 + 1.64
2.69 + 0.18
NCDS
1.37 0.40
0.98 + 0.15
2.28 + 0.24
1.18 + 0.16
Table 1: Average dynamic time warping distances (DTWD) between different approaches. NCDS exhibits comparable performance low dimensions, while clearly being superior in high dimensions.
8

---

Published as a conference paper at ICLR 2024
Perturbations
Figure 8: Robot and obstacle avoidance. Left two panels show rohot experiments with orange and red paths illustrating unperturbed and perturbed motion. The first panel shows learned vector field in R with constant orientation. In the second panel, the vector field expands TR.3 S3))), aligning the end-effector's orientation with the motion direction. Superimposed robot images depict frames of executed motion. The right two panels display how an obstacle locally reshapes the learned vector field using the modulation matrix. In the third panel, gray contours represent the learned vector field. black trajectories depict the demonstrations, and orange/green trajectories are path integrals starting from both initial points of the demonstrations and plot cornrr.. Magenta and red circles indicate the initial points of the path integrals and the obstacle, correspondingly. The last panel shows how the robot successfully avoids the obstacle (orange cylinder).
Firstly, we demonstrated 7 sinusoidal trajectories while keeping the end-effector orientation constant, which means that the robot end-effector dynamics only evolved in R As Fig. 8-a shows, the robot was able to successfully reproduce the demonstrated dynamics by following the learned vector field encoded by NCDS. Importantly, we also tested the extrapolation capabilities of our approach by introducing physical perturbations to the robot end-effector, under which the robot satisfactoriyy adapted and completed the task (see the video attached to the supplementary material).
The next experiment was aimed testing our extension to handle orientation dynamics by learning dynamics evolving in R33 x 80)3,, i.e. full-pose end-effector movements. We collected several V-shape trajectories table surface in which the robot end-effector always faces the direction of the movement, as shown in Fig. 8-b. To evaluate the performance of this setup on a more complex dataset, generated synthetic LASA dataset in R x 80(3). The results show that NCDS adeptly learns a contractive dynamical system from this dataset (see Appendix A. I). Another set of experiment was aimed at validating our approach on joint-space dynamic motions in a simulation environment (Corke & Haviland, 2021). These results show that among the baseline methods, NCDS outperforms others in generating motion within the joint space. (see Appendix 6.6.6). Furthermore, we have further extended the joint space experiments to learning human motion skills. The findings indicate that NCDS successfully learns a contractive dynamical system from complex demonstrations in high-dimensional spaces. (see Appendix A.6.2).
DISCUSSION
In this paper, we proposed the first mechanism for designing neural networks that are contractive. This allows us to build flexible class of models, called neural contractive dynamical systems (NCDS), that can learn dynamical systems from demonstrations while retaining contractive stability guarantees. Since contraction is built into the neural network architecture, the resulting model can be fitted directly to demonstration data using standard unconstrained optimization. To cope with the complexities of high sigh-dimensional dynamical systems, we further developed a variant of VAEs that ensures latent contrastive dynamical systems decode to contractive systems. Empirically, we found that the latent representation tends to behave simpler than the original demonstrations, which further increases the robustness of latent NCDS. We have further extended NCDS to the Lie group consisting of rotation matrices in order to model real-world robotic motions. Finally, we showed that our approach can be adapted to avoid obstacles while remaining contractive using the modulation techniques from Huber et al. (2022).
find that NCDS is sensitive to the choice of numerical integration scheme when solving Eq. 4. Specifically, it is important to use integration schemes with adaptive step sizing for the method to be reliable. In all experiments. we used the standard DOPRI method. One concern with adaptive step sizing is that computation time is not fixed, which may negatively influence real-time implementations. On related note. the need for numerical integration does increase the running time associated with our neural architecture compared to standard feedforward networks (see Appendix A.6.3). We, however, consider that a small price for stability guarantees.
Limitations. Empirically,

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
This work was supported by a research grant (42062) from VILLUM FONDEN. This project received funding from the European Research Council (ERC) under the European Union's 's Horizon 2020 research and innovation programme (grant agreement 757360). The work was partly funded by the Novo Nordisk Foundation through the Center for Basic Machine Learning Research in Life Science (NNF200C00626066.
REFERENCES
H. Beik-Mone Mohammadi, S. Hauberg, G. Arvanitidis, G. Neumann, and L. Rozo. Learning Riemannian manifolds for geodesic motion skills. In Robotics: Science and Systems (R:SS), 2021, URL https:
// wrrooticcoonfeeeet
org. http:///porraaaappprrrraaa /02//i/e//////eeee index...... x html.
A. Billard, S. Calinon, and Dillmann. Learning from humans. B. Siciliano and O. Khatib (eds.), Handbook of Robotics, chapter 74, pp. 1995-2014. Springer, Secaucus, NJ, USA, 2016. URL https:/ //oii. www.ro.gg...o.o.... 1007 978 3-319- --19-325522 1015522177,,, 2nd Edition.
Caroline Blocher, Matteo Saveriano, and Dongheui Lee. Learning stable dynamical systems using contraction theory. IEEE International Conference on Ubiquitous Robots and Ambient Intelli- gence (URAI), pp. 124-129, 2017. URL https: httpe/eeeeeeeeeeeeeeeee//ee ieee. org document/
7992901.
Johann Brehmer and Kyle Cranmer. Flows for simultaneous manifold learning and den- sity estimation. In Neural Information Processing Systems (NeurIPS, pp. 442-453, 2020. URL https: httts:://proceedinggs file/051 92834
www.riipscccppppe....e fwwefflle/pppe (aless/pper/22220 .pdf.
Nutan Chen, Alexej Klushyn, Alexandros Paraschos, Djalel Benbouzid, and Patrick Van der Smagt. Active learning based on data uncertainty and model sensitivity. In JEEE/RS) International Conference Intelligent Robots and Systems (IROS), pp. 1547-1554, 2018a. URL http:
///oi. 1ttp//001199//ROS.. 2018 2018.859352..
Ricky T. Q. Chen. torchdiffeq. 2018. torchdiffeq.
URL https://github... cowh www.ccooo/tggtteeeeeeoee chen/
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In Neural Information Processing Systems (NeurIPS), 2018b.
URL https://proceetings weweiiiinnnn.nuuuisseeeeeeee.eee 69386f 63388
neurips wwtpp//ppp/p./////e/eep ffie/ pdf.
Peter Corke and Jesse Haviland. Not your grandmother's toolbox the robotics toolbox reinvented for Python. In IEEE International Conference on Roborics and Automation (ICRA), pp. 11357-11363, 2021. URL https::////....... /11. 1109/ S..00/ICCA48806.00 2021 9561366.
Charles Dawson, Sicun Gao, and Chuchu Fan. Safe control with learned certificates: A survey of neural Lyapunov, barrier, and contraction methods for robotics and control. IEEE Transactions on Robotics (T-RO), 39(3):149-11767, 2023. URL https: http:/dii..oggg///// 100...1/TTOT/020222 3232542.
Luca Falorsi, Pim de Haan, Tim R. Davidson, and Patrick Forré. Reparameterizing distributions on Lie groups. International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 3244- 3253, 2019. URL https: httpp::/proceedi::: /ptp//eeeeieegs.ee//e/pppo mllr presss http://///f/lla// falorsils htmi.
Lukas Huber, Aude Billard, and Jean-Jacques Slotine. Avoidance convex and concave obstacles with convergence ensured through contraction. JEEE Robotics and Automation Letters (RA-L), 42):1462-1469, 2019, doi:
Lukas Huber, Jean-Jacques Slotine, and Aude Billard. Avoiding dense and dynamic obstacles in enclosed spaces: Application to moving in crowds. IEEE Transactions on Robotics. 38(5): 31 3113-3132, 2022. URL. https:::: wtthttps:/iiiiie //te://eeeeeeeeeeeeeeeeeee ieee. wwto (0ccummmtt/975597.
10

---

Published as a conference paper at ICLR 2024
Chia-Man Hung, Shaohong Zhong, Walter Goodwin, Oiwi Parker Jones, Martin Engelcke, loannis Havoutis, and Ingmar Posner. Reaching through latent space: From joint statistics to path planning in manipulation. JEEE Robotics and Automation Letters (RA-L), 3334533333341,, 2022. URL
http:/ /doii .rgg /t.g//0.11 1109 RTR10///RR...... 2002. 3122777
Jerome Jouffroy and Thor I. Fossen. A tutorial on incremental stability analysis using contraction theory. Modeling, Identification and Control, 2010. ISSN 0332-7353. URL https: /doi. org/10 4003.4133/mmm 2010.3.2...
S. Mohammad Khansari-Zadeh and Aude Billard. Learning stable nonlinear dynamical systems with Gaussian mixture models. IEEE Transactions on Robotics, 27(5):943-957, 2011. doi: 10.1 (100//TR..001.2155559
S. Mohammad Khansari -Zaeeh and Aude Billard. Learning control Lyapunov function to ensure stability of dynamical system-base robot reaching motions. Roboties and Autonomous Systems (RAS), 62(6)7752-765, 2014. doi:
Diederik P. Kingma and Max Welling, Auto-encoding variational Bayes. International Conference on Learning Representations (ICLR), 2014. URL https: /openrevieww net ftrm/ooormmff 33X9Ed2- 9Fy2d.
Leo Kozachkov, Michaela M Ennis, and Jean-Jacques Slotine. RNNs of RNNs: Recursive construc- tion of stable assemblies of recurrent neural networks. Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), Advances in Neural Information Processing Systems, 2022.
URL https:: /wwpeeeeeiiiieeeeeeeeee net forum? wwdddd2d3B3BBB. 38geEU..
Leo Kozachkov, Patrick Wensing, and Jean-Jacques Slotine. Generalization as dynamical robustness- the role of riemannian contraction supervised learning. Transactions on Machine Learn- ing Research, 2023. ISSN 2835-8856. URL https: http:///ppeeeeee..eemm
/ net forum?id-
f56p5mcefw.
A. Lemme, Y. Meirovitch, M. Khansari-Zadeh, T. Flash, A. Billard, and J. J. Steil. Open-source benchmarking for learned reaching motion generation in robotics. Paladyn, Journal of Behavioral Robotics, 6(1), 2015. URL http: ///////oii// org/10 1515 1511/pjbr--//e/// 2015-0002
Winfried Lohmiller and Jean-Jacques E. Slotine. On contraction analysis for non-linear systems. Automatica, 34(6):683-68) 696, 1998. ISSN 0005-1098. URL http:////oo..o//1111101110 10005-1098 (98) 00019
Jonathan Lorraine and Safwan Hossain. JacNet: Learning functions with structured Jaco- bians. INNF Workshop at the International Conference on Machine Learning (ICML), 2019. URL https: http//iveerttererrrereeeeee
hit.pp.iihhbbb... io/ INNE 2019/ acceptee papers/pdfs INNF _01 __paper__l__ pdf..
Ian R. Manchester and Jean-Jacques F. Slotine. Control contraction metrics: Convex and intrinsic criteria for nonlinear feedback design. IEEE Transactions on Automatic Control, 62(6))3046-3053 2017. URL https:://// http://i..ogg///// 111111111TTT 2017. 2668880.
Christian Mandery, Ömer Terlemez, Martin Do, Nikolaus Vahrenkamp, and Tamim Asfour. Unifying representations and large-scale whole-body motion databases for studying human motion. IEEE Transactions on Robotics, 32(4))796-809 2016. URL
pocument/750661
Klaus Neumann and Jochen J. Steil. Learning robot motions with stable dynamical systems under diffeomorphic transformations. Robotics and Autonomous Systems (RAS), 70:1-15, 2015. URL robot. 2015. 04. 606.
https: //oii /t.gg//...111. httt://...../.ooooooooooo.......
Alexander Luke Ian Norcliffe, Cristian Bodnar, Ben Day. Nikola Simidjievski, and Pietro Lio. On second order behaviour in augmented neural ODEs: A short summary. In Workshop the Symbiosis of Deep Learning and Differential Equations at NeurIPS, 2021. URL https: //openrevi ew net tttee//foommmmx.xppmttttt.

---

Published as a conference paper at ICLR 2024
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vitito Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library, 2019. URL htt.::tttps......t. /ttttts::/pprroeeeeee ngs neuripi wttppp.....pppppeee cc/ppppr/2019/ pdf.
Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, Jinkyoo Park, and Stefano Ermon. Torchdyn: Implicit models and neural numerical methods in pytorch. URL https:
Muhammad Asif Rana, Anqi Li, Dieter Fox, Byron Boots, Fabio Ramos, and Nathan Ratliff. Eu- clideanizing flows: Diffeomorphic reduction for learning stable dynamical systems. In Con- ference on Learning for Dynamics and Control (L4DC), pp. 630-639, 2020. URL https:
rana20a....mm //proceedings.mi www.eiiinnnnmm....ee.e..eee waa....pees//v1100 /atpes//a/ll////aa2eeee html.
Harish Chaandar Ravichandar and Ashwin Dani. Learning position and orientation dynamics from demonstrations via contraction analysis. Autonomous Robots, 43(4):897-912 2019. ISSN 0929-
5593. URL www.hhtpp/..oo...... //do: org/
8- 9758-x.
Harish Chaandar Ravichandar, Iman Salehi, and Ashwin Dani. Learning partially contracting dynamical systems from demonstrations. Conference on Robot Learning (CoRL), pp. 369-378, r 2017. URL https::/proceeed//e wtt.:sssml......./.em...mm www.ppressssss.
r1 7 html.
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic backpropagation and approximate inference in deep generative models. In International Conference on Machine Learning (ICML), pp. 1286, 2014.
Stefan Schaal, Auke Ijspeert, and Aude Billard. Computational approaches to motor learning by imitation. Phil. Trans. R. Soc. Lond. B, 358:537-547, 2003. URL http:://ei.org///
4098/rstbb2002.122
Malcolm David Shuster. A survey of attitude representation. Journal of The Astronautical Sciences, 41:439-517, 1993.
Vikas Sindhwani. Stephen Tu, and Mohi Khansari. Learning contracting vector fields for stable imi- tation learning. arXiv, 1804.04878, 2018. URL httptttt:/// htt.tttss///xxivve http://aaxvvv......eee /1804. 04878.
Joan Solà, Jérémie Deray, and Dinesh Atchuthan. A micro lie theory for state estimation in robotics. CoRR, abs/1812.00137,, 2018. URL http:::///xxxvv.r///aaaaa htt.//112..0 01537.
Dawei Sun, Susmit Jha, and Chuchu Fan. Learning certified control using contraction metric.
In Conference Robot Learning (CoRL), 2020. URL 12569.
Esteban G Tabak and Cristina Turner. A family of nonparametric density estimation algorithms. Communications on Pure and Applied Mathematics, 66(2):145-164, 2013.
Hiroyasu Tsukamoto and Soon Jo Chung. Neural contraction metrics for robust estimation and control: A convex optimization approach. IEEE Control Systems Letters, 5:21 2021. ISSN 1475. URL /tL.tttppe//ea....o////000/////. 2tw.:099LLLSYYS22...3...00 2020.3001- 64
Hiroyasu Tsukamoto, Soon-Jo Chung, Jcan-Jacques Slotine, and Chuchu Fan. A theoretical overview of neural contraction metrics for learning-based control with guaranteed stability. In JEEE Conference on Decision and Control (CDC), pp. 2949-2954, 2021a. URL https: / //ieeexplore. wwwpeeewwwor..eeeee....e org/ www.ment
Hiroyasu Tsukamoto, Soon-Jo Chung. and Jean-Jaques Slotine. Contraction theory for nonlinear stability analysis and suuming-basee control: A tutorial overview. Annual Reviews in Control, 52:135-169, 2021b. ISSN 1367-5788. URL https: http:::::///././///////o//////// 2021 10.000 001.
1016 arcontrol.
12

---

Published as a conference paper at ICLR 2024
J. Urain, M. Ginesi, D. Tateo, and J. Peters. Imitationflow: Learning deep stable stochastic dynamic systems by normalizing flows. In International Conference on Intelligent Robots and Systems (IROS), pp. 5231-5237. IEEE, 2020. URL http://// papercept. . empp TROS/ / fi ess 1.400.2.2.. pdf.
Julen Urain, Davide Tateo, and Jan Peters. Learning stable vector fields on lie groups. IEEE Robotics and Automation Letters, (0l:112569-157600 2022. doi: 10.11091RA..222.221900
URL htttps://ieeexploeeeeeeee ieee http///ccmmm////955500000
Keyulu Xu, Mozhi Zhang, Jingling Li, Simon Shaolei Du, Ken-Ichi Kawarabayashi, and Ste- fanie Jegelka. How neural networks extrapolate: From feedforward to graph neural net- works. In International Conference Learning Representations (ICLR), 2021. URL https: / ///ppnreview..// net wtte//oormmmooooeoooooooomo -i- -HUHHCCCCLLLCC
Jiechao Zhang, Hadi Beik-Mohammadi, and Leonel Rozo. Learning Riemannian stable dynamical systems via diffeomorphisms. In Conference on Robot Learning (CoRL), 2022a. URL https:
Yu Zhang, Long Cheng. Houcheng Li, and Ran Cao. Learning accurate and stable point-to-poing motions: A dynamic system approach. IEEE Robotics and Automation Letters (RA-L), 7(2): 1510-1517, 2022b. URL https: / http:/i.oggg//////////oooo www.oo.g//00oo 1tt///1.111111 LRA.
13

---

Published as a conference paper at ICLR 2024
A APPENDIX
A.1 NEGATIVE DEFINITENESS OF SYMMETRIC MATRIX
Let A be real n. symmetric matrix, and let As (A+AT) denote its symmetric part. Suppose that As negative definite, which means that for all nonzero vectors e E R", we have Ag Now let w be an eigenvector of A with corresponding eigenvalue X, so that Au = Aw. Then we have up Aur w Aw = Xuo w = A/w/ Taking the transpose of the equation Aur = Xun, we Xur Multiplying on the left by and using the fact that A is symmetric, we obtain w Aw = Xw" w. Therefore, we have w Aw = w Agw + w (A As) w S w "Asw < O. Since w Aw 0 for any nonzero eigenvector of A. we conclude that A is negative definite. Thus, the negative definiteness of the symmetric part As implies the negative definiteness A.
have w A
A.2 LIE OPERATORS OF (3(
To map back and forth between the Lie group 83(3) and its associated Lie algebra so(3), employ the logarithmic Log 80(3) + so(3), and exponential Exp so(3) - 80(3) maps, which are defined as follows,
Exp([r)x) =I+ sin(s)
1 -coo(()
[r|x
[r]
Ç²
R R
Log( R
(10)
2 (i)
where ( = arccos (T((R)- Moreover, R represents the rotation matrix, and [r]x denotes the 2
skew-symmetric matrix associated with the coefficient vector r.
A.3 OBSTACLE AVOIDANCE VIA MATRIX MODULATION
We here detail the used matrix modulation technique, due to Huber et al. (2022), which we use for obstacle avoidance. This approach locally reshapes the learned vector field in the proximity of obstacles, while preserving contraction guarantees. To avoid obstacles effectively, it is imperative to know both the position and geometry of the obstacle. This makes obstacle avoidance on the VAE latent space particularly difficult we need to map the obstacle location and geometry to the latent space 2. Hence, we directly apply the modulation matrix to the data space X the decoded dynamical system.
Formally, given the modulation matrix G, we can reshape the data space vector field to dynamically avoid an obstacle as follows,
= =(x)) (z) with G(x) D())
(11)
where E(x) and ((x) are the basis and diagonal eigenvalue matrices computed as, and ((x) (x=))) (x)))
E(x) = (r()) el (x) (a))..... (x) (x)],
(12)
where r(z) -- =, a reference direction computed w.r.t. a reference point æ, on the obstacle, and the tangent vectors e; form an orthonormal basis to the gradient of the distance function [(x) (see (Huber et al., 2022) for its full derivation). Moreover, the components of the matrix D are
defined as (xx) = (11- [(x] Xx(x
( (Ttx)) =1+
where p E4 a reactivity factor. Note that the matrix D modulates the dynamics along the directions of the basis defined by the set of vectors r(w) and e(x). As stated in (Huber al., 2022), the function I(-) monotonically increases w.r.t the distance from the obstacle's reference point x,, and it is, at least, a C² function. Importantly, the modulated dynamical system v = G(z) Jue (z) still guarantees contractive stability, which can be proved by following the same proof provided in (Huber et al., 2019, App. B.5).
A.4 IMPLEMENTATION DETAILS
Both the Variational Autoencoder (VAE) and the Jacobian network Jo were implemented using the PyTorch framework (Paszke al., 2019). For the VAE, employed injective generator based
14

---

Published as a conference paper at ICLR 2024
Algorithm 1: Neural Contractive Dynamical Systems (NCDS): Training in task space
Data: Demonstrations: Fn = {( R,),. E [1, E [w,
Result: Learned contractive dynamical system Tn.t Log( Rm,,): // Obtaining skew-symmetric coefficients via the Lie algebra
Pn.t [. n ::; // Create new state vector
argming_ELBB :n..): //
Train the VAE
Fn,t - //4 :); // Encode all poses using the trained VAE
=ntt = 1n,t-112m.t // Estimate latent velocities via finite differences
of
; a (0"; Prta); // Train the Jacobian network
Algorithm 2: Neural Contractive Dynamical Systems (NCDS): Robot Control Scheme
Data: Current state of the robot end-effector at time [a,, R,]
Result: Velocity of the end-effector current time step
T: Mog(R.) // Obtaining akew-symmetric coefficients
P: [zt; re]; Create new state vector
Zt Me (p:); // Compute the latent state
&, f(zz); Compute the latent velocity
J ug (zt) = i // Compute the Jacobian of the decoder
02,
d. = Jue (z)) &a; // Compute input space velocity
on M-flows (Brehmer & Cranmer, 2020), specifically using rational-quadratic neural spline flows. These flows were structured with three coupling layers. Within each coupling transform, half of the input values underwent elementwise transformation using a monotonic rational -uaadratc spline. The parameters of these splines were determined through a residual network comprising two residual blocks, with each block consisting single hidden layer containing 30 nodes. Throughout the network, we employed Tanh activations and did not incorporate batch normalization or dropout techniques. The rational -quarral-quadratic splines were constructed with ten bins, evenly distributed over the range of (-10, 10). The Jacobian network was implemented using neural network architecture consisting of two hidden layers, with each layer containing 500 nodes. The output of this network reshaped form square matrix. To facilitate the integration process, we used an off-the-shell numerical integrator called 'ddint' from the torchdi ffeq Python package (Chen, 2018 8). This package provides efficient implementations of various numerical integration methods. In experiments, we specifically utilized the Runge-Kutta and dopri5 integration methods, which are well-known and widely used for solving ordinary differential equations.
Loss functions:
The VAE is trained using the standard evidence lower bound (ELBO):
ELBO = E (x(()) [log(p((<>))) z)) KT. (=<())))\[(<>>)
(13)
where KL denotes the Kullback-Leiblrr divergence. Thereafter, the Jacobian network is trained according to the loss function Lise. which measures the error between the observed and approximated next state,
L
|z+1 (z, + 2
(14)
where 2,, 2,+1, and & represent the current and next observed latent states, along with the calculated latent velocity.
In this section, we elaborate on the training and robot control schemes for NCDS. The steps for training the Variational Autoencoder (VAE) and Jacobian network are outlined in Algo- rithm 1. Furthermore, the steps for employing NCDS to control a robot are detailed in Algorithm 2. As the algorithms show, the traning of the latent NCDS is not fully end-to-end. the latter case, we first train the VAE (end-to-end), and then train the latent NCDS using the encoded data.
Algorithms:
15

---

Published as a conference paper at ICLR 2024
0.5 1.
53
A21 31
1.0
0.5
PVLL
2205
0.
0.0
0.0
0.5
-0.5
10.0
0.0
0.0
RV3 20
-4.3
-0.5
AR
00.00 050 10
1.0
0000 0.5
20
40
60
100
120
140
20
40
60
80
100
120
140
20
40
60
80
100
120
140
Step
Stop
Step
Figure 11: Time series of demonstrations and path integrals R x 80(3)
A.5
AN INJECTIVE DECODER ARCHITECTURE OVER THE X-BALL B,
To ensure that a decoder neural network has outputs over the и-ball (Sec. 3.3), introduce simple layer. Let f IRd + RD be an injective neural network, where injectivity only considered over the image of f. If we add Tanel-laeer then the output of the resulting network the [-1, 1j box, i.e. tanh(f(=)) Rdd We can further introduce the function,
2
b(x)
|
(15)
This function smoothly (and invertible) deforms the 1] box into the unit ball (Fig. 9). Then, the function,
f(z) =b(tanh(f()))
(16)
Figure 9: An illustration of the function b(x=) in two dimensions.
injective and has the signature f Rd (D)
A.6 ADDITIONAL EXPERIMENTS
A.6.1 LEARNING 3 X SO(3)
The orientation trajectories the robot experiments tend to show relatively simple dynamics. Consequently, the learning of these tra- jectories may not fully showcase the true potential of the approach. As result, conducted additional experiments demonstrate the capacity of NCDS to encode full-pose (i.e., position-oriennttion) dynamic motions, which present more challenging and compre- hensive evaluation scenario. To construct a challenging dataset on S3(3). we projected the LASA datasets onto 3-sphere, thereby gen- erating synthetic quaternion data. Then we transformed this dataset rotation matrices on (3)).. Later, we applied the Log map project these points onto the Lie algebra. To construct the position data in 3 we simply project the 2D LASA dataset R&3 by con- catenating zero to the vector. Consequently, we concatenate this vector with the vector on the Lie algebra representing the orientation. This will produce a 6 dimensional state vector in R To increase the level of complexity, we employed two distinct LASA datasets for the position and orientation dimensions. Figure 10, depicts the latent dynamical system, the black dots represent the demonstration points, yellow path integrals start at the initial point of the demonstration, and green path integrals start in
Figure 10: The background de- picts the contours of the latent vector field, green and yellow curves depict the path integrals, and black dots represent demon- stration in latent space.
16

---

Published as a conference paper at ICLR 2024
Figure 12: Left: Latent path integrals with NCDS trained in joint space R44 Right: I Latent path integrals with NCDS trained on full human motion space (R x R$ x 80(3)). The background depicts the contours of the latent vector field, green and yellow curves depict the path integrals, and black dots represent demonstration in latent space.
random points around the demonstrations. Figure 11 shows a time series of input space trajectory in R and 83)... Moreover, the black curves show the demonstration data and other trajectories show the path integrals in each dimension. The results show that NCDS is able to learn complex non-linear contractive dynamical systems in R x (3)...
A.6.2 LEARNING HUMAN MOTIONS
To assess the model's performance in more complex environment, we conducted recent experiment using the KIT Whole-Body Human Motion Database (Mandery et al., 2016). This dataset captures subject-specifet motions, which are standardized based on the subject's 's height and weight using reference kinematics and dynamics model known as the master motor map (MMM).
Motion in Joint space R+: Focusing on the Tennis forehand motion, we selected this specific skill from the dataset, where each point in the motion trajectories is defined in 44-dimensional joint space. Figure 12-left displays the latent contractive vector field generated by NCDS. The background depicts the contours of the latent vector field, and green and yellow curves represent latent path integrals when the initial point is respectively distant from and coincident with the initial point of the demonstration. Figure 13 illustrates how the human dummy replicates each of the generated motions. The upper row of this plot depicts the progressive evolution of the demonstrated motion from left to right. Meanwhile, the middle and bottom rows illustrate the motions generated by NCDS when the initial point is respectively distant from and coincident with the initial point of the demonstration.
Motion in on full human motion space R44 R 80()): The previous experiment has primarily concentrated reconstructing motions within the joint space alone, which may be insufficient for some particular human motions. The data recorded from the human also includes the base (hip) link pose of the human dummy, which is essential for complete representation of the motion. To address this, we introduced the Kick motion skill into our experiments, where each point of the motion trajectory is defined in 44-dimensional joint space, along with the position and orientation data for the base link (hip link). This approach extends the representation the motion to comprehensive 50-dimensional input space.
Figure 12-right shows the latent path integrals generated by NCDS. Figure 14 shows the human dummy mimicking each generated motion. The upper row of this plot depicts the progressive evolution of the demonstration motion from left to right. Meanwhile, the middle and bottom rows illustrate the motions generated by NCDS when the initial point is respectively distant from and coincident with the initial point of the demonstration. This result demonstrates the inherent ability of Latent NCDS to compress a 50-dimensional nonlinear motion into 2-dimensional space.
A.6.3 QUANTITATIVE RESULTS
Regularization: It is imperative to clarify that the regularization term in Equation 2, known as epsilon plays crucial role in defining the upper bound for the eigenvalues of the Jacobian J (z).
17

---

Published as a conference paper at ICLR 2024
20
40
60
80
100
Figure 13: Generated Human Motion with NCDS trained in joint space R The upper row depicts the progressive evolution of the demonstration motion from left to right. Meanwhile. the middle and bottom rows illustrate the motions generated by NCDS when the initial point is respectively distant from and coincident with the initial point of the demonstration.
which inherently governs the contraction rate of the system. Figure 15 illustrates the effect of on the accuracy the system. Adjusting this value, particularly through increments, the potential to negatively affect the network's 's expressiveness, therefore compromising the reconstruction process. The observed negative effect can be attributed the fact that this regularization is applied to the entire Jacobian matrix and not exclusively to the elements responsible for preserving the negative definiteness of the Jacobian. As result, it might inadvertently penalize essential components of the the model's expressivity. A more targeted approach
Jacobian, leading to adverse consequences that focuses solely on the relevant elements could potentially mitigate these undesired effects and enhance the regularization's 's overall effectiveness.
To illustrate the influence of having a negative definite Jacobian in NCDS, contrast with system identical to NCDS except that remove the definiteness constraint on the Jacobian network. Figure 16 shows the path integrals generated by both contractive (Figure 16-a) and unconstrained (Figure 16-b) dynamical systems. As anticipated, the dynamics generated by the unconstrained system lack stable behavior. However, the vector field aligns with the data trends in the data support regions. Furthermore, we performed similar experiments with two different baseline approaches: an MLP and Neural Ordinary Differential Equation (NeuralODE) network, to highlight the effect of the contraction constraint on the behavior of the dynamical system. The MLP baseline uses neural network with 2 hidden layers each with 100 neurons with tanh activation function. The NeuralODE baseline is implemented based on (Poli et al.), with 2 hidden layers each with 100 neurons. The model is then integrated into a NeuralODE framework. The NeuralODE is configured with the 'adjoint'' sensitivity method and utilizes the 'dopri5' solver for both the forward and adjoint passes. is important to mention that these baselines directly reproduce the velocity according to = fz) These models were trained with the same cost function as our NCDS's Jacobian network loss defined in Equation 14. Both networks trained for 1000 epochs
Unconstrained dynamical system:
18

---

Published as a conference paper at ICLR 2024
25
100
Figure 14: Generated Human Motion with NCDS trained on full human motion space (R 44 x R x 80(3): The
upper row depicts the progressive evolution of the demonstration motion from left to right. Mcanwhile, the middle and bottom rows illustrate the motions generated by NCDS when the initial point is respectively distant from and coincident with the initial point of the demonstration.
() 10
(b) 10
10
(d)
Figure 15: Path integrals generated by NCDS trained using different regularization term e.
with ADAM optimizer. Figure 16-c and 16-d shows that both models effectively model the observed dynamics (vector field); however, anticipated, contraction stability is not achieved.
Furthermore, we tested the integrator in a toy example constructed using a Sine wave under different conditions. Figure 17 shows the integrator under different conditions. Figure 17-a and -b displays cases when the absolute and relative tolerance parameters the integrator need to be at least as 10 Figure 17-c shows that the number of time steps to solve the integration does not make any difference. Furthermore, Figure 17-d displays different integration methods and their data reconstruction accuracy.
Integration:
19

---

Published as a conference paper at ICLR 2024
(a) Contractive
(b) Unconstrained
MLP
(d) Noural
Figure 16: Path integrals generated under the Neural Contractive Dynamical Systems (NCDS) setting, along with baseline comparisons using Multilayer Perceptron (MI.)) and Neural Ordinary Differential Equation (NeuralODE) models.
4
322
(2) Comparison umong different absolute tolerance values
(b) Comparison amoog different relative tolerunce values
t0 shos
Steps
Aestluu
Comparison among different integration time steps
Comparison among different integration methods
Figure 17: Evaluation of the integrator
Activation function: The choice of activation function certainly affects the generalization in the dynamical system (where there is no data). To show this phenomena, we ablated few common activation functions. For this experimental setup, we employed a feedforward neural network comprising two hidden layers, each consisting of 100 units. As shown in Figure 18, both the Tanh and Softplus activation functions display superior performance, coupled with satisfactory generalization capabilities. This indicates that the contour of the vector field aligns with the overall demonstration behavior outside of the data support. Conversely, opting for the Sigmoid activation function might yield commendable generalization beyond the confines of the data support. However, it is essential to acknowledge that this choice could potentially encounter challenges when it comes to reaching and stopping the target.
Dimensionality and execution time: Table 2 compares 5 different experimental setups. First, Table 2-a reports the average execution time under the setup where the dimensionality reduction (i.e., no VAE) is discarded and the dynamical system is learned directly in the input space. As the results show, increasing the dimensionality increase in execution time for single integration step. This empirical relationship reaffirms the inherent computational intensity tied Jacobian-based computations. approach leverages the dimensionality reduction provided by custom Variational Autoencoder (VAE) that takes advantage an injective generator (as decoder) and its inverse (as encoder). As showcased in Table 2-b, the VAE distinctly mitigates the computational burden. Comparing Table 2-a and Table 2-b, the considerable advantages of integrating the contractive dynamical system with the VAE into the full pipeline become evident it leads to significant halving of the execution time in the 8D scenario. Additionally, the findings indicate that increasing the dimensionality from 8 to 44 results in sixfold increment in execution time. Moreover, the presented results indicate that
the input space from 2 entails remarkable 40-fold
address this issue, our
20

---

Published as a conference paper at ICLR 2024
(n)
(b) Softplus
(e) Sigmoid
Figure 18: Path integrals generated by NCDS trained using different activation functions.
Input. Dim
3D
8D
44D
Input. Dim 2D
8D
time (ms)
10.6
20.2
121.0
time (ms)
1.23
40.9
Learning the vector field in the latent space (full pipeline).
Learning the vector field directly in the input space.
Table 2: Average execution time milliseconds) a single integration step
current system may not achieve real-time operation. However, it is noteworthy that the video footage from our real-world robot experiments convincingly demonstrates the system's 's rapid query response time, sufficiently efficient control the robot arm, devoid of any operational concerns.
A.6.4 COMPARISONS
Table 3 summarizes the comparisons presented in the main papers in the contraction literature. Here we note lack consensus regarding the method or baseline for comparison the context of learning contractive dynamical systems. A significant number of approaches have overlooked the contraction properties of their baseline, opting instead to base their comparisons solely on asymptotic stability guarantecs. Following this trend, have conducted thorough and extensive comparative analysis involving Stable Estimator of Dynamical Systems (SEDS), Euclideanizing Flows (EF), and Imitation Flows (IF). We emphasize that neither of these baselines provides contraction guarantees but rather asymptotic stability guarantees. Therefore, the comparison is provided solely to evaluate the general stability and accuracy of the reconstructed motion.
To perform comparison between these methods and ours, we have provided two groups of datasets. The first group comprises multiple synthetic datasets, derived from the LASA dataset. The second category combines joint space trajectories of 7-DOF Tranka-Emike Panda robot arm.
Experimental setups:
Imitation Flow: For Imitation Flow we have used imitation flow network with the depth of 5. Each level within this depth was constructed utilizing CouplingLayer, Ran- domPermutation,, and LULinear techniques, as recommended by the authors for managing high-ding -igh--imenssonal data. The training process includes 1000 epochs, with empirically chosen learning rate of 10-3 These hyper parameters were found experimentally and the best results have been selected to be compared against.
Euclidenizing Flow: Here, have used coupling network with random Fourier features parameterization, with 10 coupling layers each with 200 hidden units, and length scale of 0.45 as the length scale for random Fourier features. The training process includes 1000 epochs, with empirically chosen learning rate of 10-4 These hyper parameters were found experimentally and the best results have been selected to be compared against.
SEDS: Here, we use 5 Gaussian functions. The training process includes 1000 epochs with MSE objective.
21

---

Published as a conference paper at ICLR 2024
A.6.5 SYNTHETIC LASA
First, evaluate the performance and accuracy of the methods 2D LASA dataset the baseline for the experiments. As depicted in Fig. 19, all methods demonstrate strong performance, effectively capturing the data trends and successfully reaching the target points.
NCDS
(b) SEDS
(c) EF
(d)
Figure 19: Path integrals generated using different approaches for the two-dimensional LASA problem. All methods provide similar reconstructions.
Examining the average time warping distance in Fig. 20, becomes evident that EF and IF have better performance in reproducing the system dynamics.
2.0
18
Megn
10
NCDS
SEDs
Figure 20: Average dynamic time warping distance between path integrals and the demonstrations for the two-dimensional LASA problem.
We extend our evaluation to 4D LASA synthetic dataset created by merging two distinct datasets. In Figure 21, we observe that while NCDS and maintain a strong level of performance, SEDS and EF exhibit comparatively inferior results.
Mean
15.
1.:
35
NCD5
Figure 21: Average dynamic time warping distance between path integrals and the demonstrations for the four-dimensional LASA problem.
Lastly, we assess the methods using an 8D LASA synthetic dataset, constructed by concatenating different datasets,
D Angler MAngles Tine MLine TSine VSines TiShape MiShapel E BR Tx8
(17)
As illustrated in Fig. 22, NCDS has the best performance in higher-dimensional spaces.
22

---

Published as a conference paper at ICLR 2024
10
(een DTVDD
SEbs
Figure 22: Average dynamic time warping distance between 5 trajectories constructed with different methods and demonstration trajectories
It is worth noting that in the 2D scenario, NCDS was directly trained in the input space, whereas in the 4D and 8D experiments, the dynamics were learned in the latent space
To make sure that we have a fair comparison we have scaled up all the methods. 1. Imitation flow: the depth the network was scaled up by factor of (from to 10), We evaluated 20 different architecture with higher and lower depth and used the one with the best performance. 2. Euclideanizing Flow: the number of coupling blocks in the BijectionNet in the original code released by the Authors was increased by the factor of from (10 to 20). We evaluated 20 different architecture with higher number of blocks and nodes and used the one with the best performance. SEDS: we increased the number of Gaussian functions from 5 to 20. We evaluated 10 different architecture with different number of Gaussian functions used the one with the best performance. 4. NCDS: The Jacobian network did need more parameters since the latent space has stayed 2D. Therefore, we just added the VAE to the architecture when moving from 2D to 8D.
A.6.6 JOINT SPACE DYNAMICS
Since using the pose information including Quaternion is not feasible due to the fact that none of the methods that compared against support non-Fuclidean spaces, have chosen to joint space trajectories of 7-DOF robot. Therefore, demonstration data that is used is 7-dimensional and represents "V" shape character table. Fig. 23 shows the comparison on the average dynamic time warping distance (DTWD) among different methods when learning 7D joint space trajectories.
Mean DTWD
NCDS
Figure 23: Average dynamic time warping distance between path integrals and the demonstrations for the 7-dimensional Joint space problem.
As shown, NCDS provides the highest reconstruction accuracy, therefore showing how our method better scales higher-dimensional settings. The video of this experiment int the supplementary material. Moreover, the video shows the recording process of the demonstrations and also 4 dif- ferent videos from simulated robots performing the path integrals in the joint space using the four benchmarked methods.
23

---

Published as a conference paper at ICLR 2024
Paper title
Abbr.
Compared against
Year
Control with Learned Certificates: A Survey Neural Lyapunov, Barrier, and Contraction methods
2022
Learning Contraction Policies from Offline Data
MPC (I.QR)) (CQT.)
2022
Neural Contraction Metries for Robust Es- NCM
CV.STEM, LQR
2021
timation and Control: Convex Optimiza- tion Approach
Learning Stabilizable Nonlincar Dynamics
CCM-R
2021
with Contraction- Based Regularization
Learuing Certified Control Using Contrac-
C3M
SoS (Sum-of-Squares programming), MPC,
2020
tion Metric
RL (PPO)
Learning position and orientation dynamics
CDSP
SEDS. CLF-DM (Control Lyapunov
2019
from demonstrations via contraction analy-
Function-based Dynamic Movements),
Tau-SEDS, NIVF (ncurally imprinted
vector fields)
Leuminng Stable Dynamical Systems using
C-GMR
SEDS
2017
Contraction Theory
Learning Contracting Voctor Fields For Sta- CVF
DMP, SEDS and CLE-DM
2017
ble Imitation Learning
Leamning Partially Contracting Dynamical
CDSP
(only
DMP. CLF-DM
2017
Systems from Demonstrations
position)
Table 3:
The present state-of-the-art literature pertaining to the learning of contractive dynamical systems.
24


Reviews:
[{"review_id": "KoGRKP6tjp", "reviewer": "Reviewer_4cHX", "summary": "This paper introduces a method to learn contractive dynamical systems. One method of constructing contractive systems has been to apply a diffeomorphism to provide a change of coordinates to a known contractive system. This paper proposes to use a VAE instead of a diffeomorphism, and constructs the decoder such that it is injective. As such, by enforcing that the latent dynamics is contractive, the dynamics in the data space will also be contractive. The paper also extends the method to Lie groups to account for end-effector orientation.\n\nThe paper is generally well-written, the method appears sound, and the motivations are clear. My main concerns are as follows:\n\n1. Theorem 1 states that the contractivity is perserved under a diffeomorphism. There is a lack of analysis and formal guarentees around the assumption that an injective function acting on a lower dimensional contractive system produces a higher-D contractive system. Just as a thought experiment, what happens to coordinate points in the higher-D data space where there does not exist a coordinate in the lower-D latent space? \n\n2. How is the collision-avoidance on the entire manipulator handled? The learned dynamical system seems to model the end-effector pose, but collision-avoidance should be handled across the body of the robot. One approach to handle this is to pull the dynamical system to the C-space of the robot and define body-points for the collision-avoidance, as done in (Zhi 2022, L4DC). This is a relevant reference and should be reviewed, as it also takes a diffeomorphic learning approach.\n\nOverall, I believe this is a neat idea, but more clarity around the theoretical insights is needed. I'm happy to raise my score when my concerns have been address.", "questions": "See above.", "limitations": "", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "See above.", "weaknesses": "See above.", "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": "dEydy6FGQB", "reviewer": "Reviewer_fKvz", "summary": "This work addresses the challenging problem of learning contractive dynamical systems using neural networks. The key idea is to utilize the fact that contractivity is invariant under diffeomorphisms. This motivates the use of autoencoders to learn contractive dynamics in a low-dimensional latent space. The proposed method includes a VAE architecture which naturally enforces the contractivity of the dynamics in the latent space. These results are then applied to several interesting applications, including obstacle avoidance, and learning dynamics in $SO(3)$.", "questions": "* This work focuses on learning contractive dynamics. However, suppose the system $\\dot{x} = f(x)$ is *not* contractive. What can you say about the solution to the optimization problem (i.e. the feasibility of eqns (4) and (5)) in this scenario?\n* In a similar vein, it seems like this work could easily be extended to (neural) controller synthesis in the latent space. Can the authors comment on this?\n* Can the authors comment on the practical effectiveness of the model in greater detail? As mentioned in the discussion section (p9), the cost of numerical integration can be extensive. Have the authors come across examples where this has been an impediment to performance?\n* Could the authors also address the concerns raised in the weaknesses section?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The paper has several strengths.\n\n* The idea of learning a low dimensional latent space embedding of the dynamics is interesting and novel, with a variety of interesting potential applications\n* The construction of the VAE ensures that the contractivity is invariant to the mapping to the latent space (and vice versa)\n* The proposed framework is extended to a variety of scenarios, including dynamics over Lie groups ($SO(3)$) and obstacle avoidance\n* The paper itself is, generally speaking, well written.", "weaknesses": "There are a few weaknesses.\n\n* It would be nice to have the invariance of the contractivity stated formally. \n* In the discussion section (page 9), it is mentioned that the choice of integration scheme can siginificantly affect the behaviour of the learned model. This requires further discussion. For instance, how significantly does the computation time affect the performance of the model? Is there a choice of integrator that doesn't require adaptive step-sizes (perhaps a symplectic integrator)?\n* It would be helpful if the training algorithm were stated explicitly in Section 3.", "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": "OizrLboDSf", "reviewer": "Reviewer_2Vh1", "summary": "### Problem Statement\nThe paper discusses the significant problem concerning the assurance of stability of data-driven controlling of robots, especially when the learned dynamics are controlled by neural networks. Stability is critical to prevent robots from executing harmful or undesirable actions. However, achieving global stability in dynamical systems learned from data proves challenging.\n\n### Main Contribution\nThe primary contribution of this paper is the introduction of a new methodology called Neural Contractive Dynamical Systems (NCDS), which can guarantee contractive stability for the dynamics it learns, in both Euclidean and SO(3) manifolds.\nThis makes NCDS a highly adaptable learning architecture offering contractive stability guarantees and obstacle avoidance capabilities. The empirical results show that their approach encodes desired dynamics more accurately compared to existing state-of-the-art methods while providing stronger stability guarantees. Through NCDS, the authors aim to bridge the gap between learning robot dynamics from demonstrations and ensuring stability, which has been a notable challenge due to the extrapolating behavior of neural network models.\n\n### Methodology\nThis methodology is designed to learn stable dynamical systems by constructing negative definite Jacobian from the output of neural network to ensure contraction, thus leading to global stability. The authors extend this method to address high-dimensional dynamical systems by developing a variant of the variational autoencoder using flow-based diffeomorphisms, which learns dynamics in a low-dimensional latent space while maintaining contractive stability post-decoding. They further extend their methodology to include contractive systems on the Lie group of rotations, catering to full-pose end-effector dynamic motions. The method can also be incorporated with a matrix-modulation technique to enable obstacle avoidance.\n\n### Experiments\nThe paper evaluates the effectiveness and scalability of Neural Contractive Dynamical Systems (NCDS) through synthetic and real-world tasks. Initial tests on 2D trajectories from the LASA dataset demonstrated NCDS's capability in capturing and replicating underlying dynamics, even in regions not covered by original data. Compared to baseline methods like Euclideanizing flow, Imitation flow, and SEDS, NCDS was the only method showing contractive behavior indicative of stability. When scaled to higher-dimensional data (LASA-4D and LASA-8D datasets), NCDS maintained a good performance, contrary to the deteriorating performance of baseline methods.\n\nFurthermore, the obstacle avoidance capability of NCDS was showcased on the LASA dataset, with successful generation of safe trajectories around obstacles. Real-world robot experiments on a 7-DoF Franka-Emika robotic manipulator underlined NCDS's effectiveness in reproducing demonstrated dynamics and adapting to physical perturbations. The experiments collectively underline NCDS's potential in managing various aspects of robotic motion learning, ensuring stability, and navigating obstacles, crucial for advancing real-world robotic applications.", "questions": "1. Implementation details\n   1. Is training end-to-end? What are the loss functions used for training?\n   2. How are the sequential data used for training? Is each trajectory used as a whole as one data point, or cut into segments as multiple data points?\n   2. Second-order?\n\n1. Multi-modality", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "### Originality and Significance\n\nThis is the first work to my knowledge that endows neural network based dynamics modeling with guaranteed contractive stability, which is, as the authors point out, important to robotics as people are trying to take advantage of the modeling capacity of neural networks. The conservative and non-diverging extrapolation that comes with the contractiveness enabled by this work could also benefit neural modeling of other dynamical systems apart from robotics, as the neural ordinary differential equations are known to have difficulty modeling dynamical systems when data is not enough to cover the state space, as the unregularized extrapolation of neural networks could easily lead to numerical instability during integration. Therefore this work is not only innovative but also has potentially significant impact to the community.\n\n### Quality\n\nThe work derives the method from contraction theory and designed practical algorithms for data-driven robotics tasks, and compared with several state of the art baselines.\n\n### Clarity\n\nThe writing of the paper is great. It excellently motivates the work and explains the method and results well with proper details.", "weaknesses": "1. Some implementation details are missing and the codes are not available. See more in the Questions section.\n\n2. Limitation of baseline method choices: The baselines compared to are all focused on asymptotic stability guarantees. While these are the most relevant methods for comparison, it would be interesting to see how imitation learning methods without any stability guarantee works on the tasks, to demonstrate the necessity of stability guarantee.", "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 (spotlight)