Learning Neural Contracting Dynamics: Extended Linearization and Global Guarantees

Authors: sean  jaffe
alexander  davydov
deniz  lapsekili
ambuj k singh
francesco  bullo
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 01:39:57

Abstract: Global stability and robustness guarantees in learned dynamical systems are essential to ensure well-behavedness of the systems in the face of uncertainty. We present Extended Linearized Contracting Dynamics (ELCD), the first neural network-based dynamical system with global contractivity guarantees in arbitrary metrics. The key feature of ELCD is a parametrization of the extended linearization of the nonlinear vector field. In its most basic form, ELCD is guaranteed to be (i) globally exponentially stable, (ii) equilibrium contracting, and (iii) globally contracting with respect to some metric. To allow for contraction with respect to more general metrics in the data space, we train diffeomorphisms between the data space and a latent space and enforce contractivity in the latent space, which ensures global contractivity in the data space. We demonstrate the performance of ELCD on the high dimensional LASA, multi-link pendulum, and Rosenbrock datasets.

Main paper: 
 introduction: 
 Due to their representation power, deep neural networks have become a popular candidate for modeling continuous-time dynamical systems of the formEQUATIONwhere f (x) is an unknown (autonomous) vector field governing a dynamical process. Beyond approximating the vector field f , it is desirable to ensure that the learned vector field is well-behaved. In many robotic tasks like grasping and navigation, a well-behaved system should always reach a fixed endpoint. Ideally, a learned system will still stably reach the desired endpoint even when pushed away from demonstration trajectories. Additionally, the learned system should robustly reach the desired endpoint in the face of uncertainty. In tasks such as manufacturing, animation, and human-robot interaction, functionality and safety require the learned system must smoothly follow a specific trajectory to its target.To enforce stability guarantees, a popular approach has been to ensure that the learned dynamics admit a unique equilibrium point and have a Lyapunov function, e.g. [30] . While popular, approaches based on Lyapunov functions typically struggle to provide general robustness guarantees since global asymptotic stability does not ensure robustness guarantees in the presence of disturbances. Indeed, input-to-state stability of global asymptotically stable dynamical systems needs to be separately established, e.g., see Chapter 5 in [21] .To ensure robustness and to allow for smooth trajectory following, there has been increased interest in learning contracting dynamics from data [6] . A dynamical system is said to be contracting if any two trajectories converge to one another exponentially quickly with respect to some metric [27] . If a learned contracting system that admits a demonstration trajectory is pushed off that trajectory, it will follow a new trajectory that exponentially converges to the demonstration trajectory. Additionally, if a system is contracting, it is exponentially incrementally input-to-state stable [39] . If the system is autonomous, it also admits a unique equilibrium, is input-to-state stable, and has two Lyapunov functions that establish exponential stability of the equilibrium [9] . Since establishing contractivity globally requires satisfying a matrix partial differential inequality everywhere, prior works have focused on establishing contractivity in the neighborhood of training data [35, 36] .Stable, but not necessarily contracting dynamics. Numerous works have aimed to learn stable dynamical systems from data including [7, 15, 22, 30, 32, 40, 45] . In [22] , the authors introduce the Stable Estimator of Dynamical Systems (SEDS), which leverages a Gaussian mixture model to approximate the dynamics and enforce asymptotic stability via constraints on learnable parameters. In [30] , the authors jointly learn the dynamics and a Lyapunov function for the system and project the dynamics onto the set of dynamics which enforce an exponentially decay of the Lyapunov function. This projection is done in closed-form and establishes global exponential convergence of the dynamics to the origin. In [40] , the authors introduce Imitation Flow where trajectories are mapped to a latent space where states evolve in time according to a stable SDE. In [32] , Euclideanizing Flows is introduced where the latent dynamics are enforced to follow natural gradient dynamics [2] . A similar approach is taken in [45] where additionally collision avoidance is considered.Existing works on learning contracting dynamics. Learning contracting vector fields from demonstrations has attracted attention due to robustness guarantees [6, 33, 35, 37, 39] . In [33] , the dynamics are defined using a Gaussian mixture model and the contraction metric is parametrized to be a symmetric matrix of polynomial functions. Convergence to an equilibrium is established using partial contraction theory [41] . In [35] , the dynamics are defined via an optimization problem over a reproducing kernel Hilbert space; the dynamics are constrained to be locally contracting around the data points, i.e., there may be points in state space where the dynamics are not contracting. In [37] and [39] , the authors study controlled dynamical systems of the form ẋ = f (x, u), where u is a control input and train neural networks to find a feedback controller such that the closed-loop dynamics are approximately contracting, i.e., contractivity is not enforced but instead lack of contractivity is penalized in the cost function.Recent work, [6] , has proposed a Neural Contractive Dynamical System (NCDS) which learns a dynamical system which is explicitly contracting to a trajectory. NCDS parametrizes contracting vector fields by learning symmetric Jacobians and performing line integrals to evaluate the underlying vector field. NCDS is computationally costly because of this integration. Additionally, the Jacobian parametrization used is overly restrictive, because not all contracting vector fields have symmetric Jacobians. Constraining the vector field to have symmetric Jacobian is equivalent to enforcing that the vector field is a negative gradient flow and contracting with respect to the identity metric [42] .For scalability to higher dimensions, NCDS leverages a latent space structure where an encoder maps the data space to a lower-dimensional space and enforces contracting dynamics in this latent space. Then a decoder "projects" the latent-space dynamics to the full data space. It is then argued that on the submanifold defined by the image of the decoder, the dynamics are contracting.In this paper, we present a novel model for learning deep dynamics with global contraction guarantees.We refer to this model as Extended Linearized Contracting Dynamics (ELCD). To the best of our knowledge, ELCD is the first model to ensure global contraction guarantees. To facilitate the development of this model, we provide a review of contracting dynamics and extended linearization of nonlinear dynamics. Leveraging extended linearization, we factorize our vector field as f (x) = A(x, x * )(x -x * ), where x * is the equilibrium of the dynamics. We enforce negative definiteness of the symmetric part of A(x, x * ) everywhere and prove (i) global exponential stability of x * , (ii) equilibrium contractivity of our dynamics to x * , and (iii) using a converse contraction theorem, contractivity of the dynamics with respect to some metric.Since negative definiteness of the symmetric part of A is not sufficient to capture all contracting dynamics, we introduce a latent space of equal dimension as the data space and learn diffeomorphisms between the data space and this latent space. The diffeomorphisms provide additional flexibility in the contraction metric and allow learning of arbitrary contracting dynamics compared to those which are solely equilibrium contracting.Our example in Section 3.4 provides theoretical justification for why the diffeomorphism and the learned contracting model must be trained jointly. In summary, if the diffeomorphism is trained first, and transforms the data to a latent space, the model class may not be expressive enough to accurately represent the true latent dynamics. If the diffeomorphism and dynamics are trained simultaneously, this limitation is overcome. This is in contrast to [6] which trains the diffeomorphism independently as a variational autoencoder and then trains the model after on the transformed data.Our model, ELCD, directly improves on NCDS in several ways. We parameterize the vector field directly instead of parametrizing its Jacobian. Doing so prevents us from needing to integrate the Jacobian to calculate the vector field and thus speeds up training and inference. ELCD is also more expressive than NCDS because it can represent vector fields with asymmetric Jacobians. Additionally, ELCD is guaranteed to be contracting to an arbitrary equilibrium point, either selected or learned, at all training steps. NCDS, in contrast, must learn the equilibrium point over the course of training. Additionally, NCDS learns an encoder and decoder for a lower-dimensional latent space and thus can only be contracting on the submanifold defined by the image of the decoder. From initial conditions that are not on this submanifold, NCDS may not exhibit contracting behavior. In contrast, since the latent space of ELCD is of the same dimension as the data space, we train diffeomorphisms that ensure global contractivity in the data space. ELCD exhibits better performance than NCDS at reproducing trajectories from the LASA [26] , n-link Pendulum, and Rosenbrock datasets. We additionally compare against the Euclideanizing Flow [32] , and Stable Deep Dynamics [30] models. 
preliminaries: 
 We consider the problem of learning a dynamical system ẋ = f (x) using a neural network that ensures that the dynamics are contracting in some metric. Going forward, we denote by Df (x) = ∂f ∂x (x) the Jacobian of f evaluated at x. To this end, we define what it means for a dynamical system to be contracting. Definition 1. A contracting dynamical system is one for which any two trajectories converge exponentially quickly. From [27] , for a continuously differentiable map f : R d → R d , the dynamical system ẋ = f (x) is contracting with rate c > 0 if there exists a continuously-differentiable matrixvalued map M : R d → R d×d and two constants a 0 , a 1 > 0 such that for all x ∈ R n , M (x) = M (x) ⊤ and a 0 I d ⪰ M (x) ⪰ a 1 I d and additionally satisfies for all x M (x)Df (x) + Df (x) ⊤ M (x) + Ṁ (x) ⪯ -2cM (x).(EQUATIONThe map M is called the contraction metric and the notation Ṁ (x) is shorthand for the n × n matrix whose (i, j)entry is Ṁ (x) ij = ∇M ij (x) ⊤ f (x).A central result in contraction theory is that any dynamical system ẋ = f (x) satisfying (2) has any two trajectories converging to one another exponentially quickly [27, 29, 39] . Specifically, there exists L ≥ 1 such that any two trajectories x 1 , x 2 of the dynamical system satisfyEQUATIONIn other words, contractivity establishes exponential incremental stability [4] .We note that any contracting dynamical system enjoys numerous useful properties including the existence of a unique exponentially stable equilibrium and exponential input-to-state stability when perturbed by a disturbance. We refer to [27] for more useful properties of contracting dynamical systems and [9] for a recent monograph on the subject.The problem under consideration is as follows: given a set of demonstrations D = {(x i , ẋi )} n i=1 consisting of a set of n state, x i ∈ R d , and velocity, ẋi ∈ R d , pairs, we aim to learn a neural network f (x i ) = ẋi that parametrizes a globally contracting dynamical system with equilibrium point x * , such that f (x * ) = 0. In essence, this task requires learning both the vector field and the contraction metric, M such that they jointly satisfy (2).Suppose we assumed that the dynamical system we aimed to learn was linear, i.e., ẋ = Ax for some matrix A ∈ R d×d and we wanted to find a contraction metric M which we postulate to be constant M (x) := M for all x. The contraction condition (2) then readsEQUATIONwhich implies that A has all eigenvalues with Re(λ(A)) ≤ -c, see Theorem 8.2 in [18] . In other words, for linear systems, contractivity is equivalent to stability.Although the condition (4) is convex in M at fixed A, the task of learning both (A, M ) simultaneously from data is not (jointly) convex. Instead, different methods must be employed such as alternating minimization. A similar argument is made in [30] in the context of stability and here we show that the same is true of contractivity.In the case of nonlinear systems, contractivity is not equivalent to asymptotic stability. Indeed, asymptotic stability requires finding a continuously differentiable Lyapunov function V : R d → R ≥0 which is positive everywhere except at the equilibrium, x * , and satisfies the decay conditionEQUATIONOne advantage of learning asymptotically stable dynamics compared to contracting ones is that the function V is scalar-valued, while the contraction metric M is matrix-valued. A disadvantage of learning asymptotically stable dynamics compared to contracting ones is that we are required to know the location of the equilibrium point beforehand and no robustness property of the learned dynamics is automatically enforced. To this end, there are works in the literature that enforce an equilibrium point at the origin and learn the dynamics and/or the Lyapunov function under this assumption [30] . In the case of contractivity, existence and uniqueness of an equilibrium point is implied by the condition (2) and it can be directly parametrized to best suit the data.3 MethodsIn this section, we describe the models that we compare to, namely EFlow [32] , SDD [30] , and NCDS [6] .Euclideanizing Flow: Starting with EFlow, they parametrize their dynamical system byEQUATIONwhere ψ : R d → R d is a diffeomorphism, Φ : R d → R is a convex potential function, and G ψ (x) = Dψ(x) ⊤ Dψ(x) is the induced Riemmanian metric. The dynamics (18) then has the interpretation of being the steepest descent for Φ • ψ on the Riemmanian manifold defined by metric G ψ . Specifically, in [32] , the convex potential function is defined to be Φ(y) = ∥y -y ⋆ ∥ 2 , where y ⋆ = ψ(x ⋆ ) and the diffeomorphism is parametrized via a coupling layer [12] . Note that contraction properties of dynamics of these form were studied in [42] . Since Ψ is chosen to be convex but not strongly convex, these dynamics are only weakly contracting, i.e., satisfy (2) with c = 0.Stable Deep Dynamics: In [30] , the authors parametrize both an unconstrained vector field f : R d → R d and a Lyapunov function V : R d → R ≥0 via neural networks. Then to enforce global exponential stability, at every point x ∈ R d , they project f (x) onto the convex set of points {u ∈ R d | ∇V (x) ⊤ u ≤ -αV (x)}. They do this projection in closed-form so that the dynamics have the representationEQUATION)Under a smart parametrization of V using input-convex neural networks [3] , the authors guarantee global exponential stability of ( 19) to the origin with rate α. Note that the dynamics ( 19) are continuous, but not necessarily differentiable. Due to this potential nonsmoothness, it is theoretically unknown whether these dynamics are contracting (since Theorem 2.8 in [16] requires at least some differentiability).Neural Contractive Dynamical Systems: In [6] , the authors enforce contractivity by directly enforcing negative definiteness of the Jacobian matrix of the vector field by parametrizingEQUATIONwhere J θ : R d → R d×d is a neural network. Under this parametrization of the Jacobian, it is straightforward to see that (2) holds with c = ϵ and M (x) = I d for all x. To recover the vector field from the Jacobian, the fundamental theorem of calculus for line integrals (which can be seen as a version of the mean-value theorem) is utilized to expressEQUATIONwhere x 0 and f (x 0 ) denote the initial condition and its velocity, respectively.To express richer dynamics, NCDS also consists of an encoder and decoder for a lower-dimensional latent space. The dynamics in the latent space are constrained to be (21) and the encoder and decoder are parametrized using a VAE. Since, generally, the latent space is of lower dimension than the data space, the resulting dynamics are only guaranteed to be contracting on the submanifold defined by the image of the decoder. 
motivation: 
 The motivation for our approach comes from the well-known mean-value theorem for vector-valued mappings, which we highlight here. Lemma 2 (Mean-value theorem (Proposition 2.4.7 in [1] )). Let f : R d → R d be continuously differentiable. Then for every x, y ∈ R d ,EQUATIONIf y = x * satisfies f (x * ) = 0, then the continuously differentiable map f admits the factorizationEQUATIONEQUATION)This factorization is referred to as extended linearization in analogy to standard linearization, i.e., for x close to x * , f (x) ≈ Df (x * )(x -x * ). We remark that when d ≥ 2, extended linearization is not unique, and for a given f , there may exist several A such that f (x) = A(x, x * )(x -x * ). In other words, (8) showcases one valid choice for A such that this factorization holds. Indeed, this nonuniqueness has been leveraged in some prior works, e.g. Section 3.1.3 in [39] , to yield less conservative contractivity conditions.Since we know that a contracting vector field admits a unique equilibrium point, x * , we restrict our attention to learning a matrix-valued mapping A : R d × R d → R d×d and ensuring that this mapping has enough structure so that the overall vector field satisfies the contraction condition (2) for some contraction metric M . This task is nontrivial since fEQUATIONwhere ∂A ∂x : R d × R d → R d×d×d is a third-order tensor-valued mapping. In what follows, we will show that a more simple condition will imply contractivity of the vector field. Namely, negative definiteness of the symmetric part of A(x, x * ) will be sufficient for contractivity of the dynamical system ẋ = A(x, x * )(x -x * ).A(x, x * )We show a simple example of our model, the Extended Linearized contracting Dynamics (ELCD). Let x ∈ R d be the state variable and f : R d → R d be a vector field with equilibrium point x * . As indicated, we parametrize our vector field by its extended linearizationEQUATIONwhere nowEQUATIONP s , P a : R d × R d → R d×d are neural networks, I d is the d-dimensional identity matrix, and α > 0 is a constant scalar. Note that the symmetric part of A(x, x * ) is negative definite sinceA(x, x * ) + A(x, x * ) ⊤ 2 = -P s (x, x * ) ⊤ P s (x, x * ) -αI d and P s (x, x * ) ⊤ P s (x, x *) is guaranteed to be positive semidefinite.We prove that a vector field parametrized this way is guaranteed to be (i) globally exponentially stable, (ii) equilibrium contracting as defined in [10] , and (iii) contracting in some metric. The key tools are partial contraction theory [41] and a converse contraction theorem. Theorem 3 (Equilibrium Contraction and Global Exponential Stability). Suppose the dynamical system ẋ = f (x) is parametrized as f (x) = A(x, x * )(x -x * ) where A(x) is as in (11) . Then any trajectory x(t) of the dynamical system satisfiesEQUATIONProof. We use the method of partial contraction as in [41] . Let x(t) be a solution to ẋ(t) = A(x(t), x ⋆ )(x(t) -x * ) with initial condition x(0) = x 0 . Then, define the time-varying virtual system ẏ(t) = A(x(t), x * )(y(t) -x * ).(13) We will establish that this virtual system is contracting in the identity metric, i.e., ( 2) is satisfied with M (x) = I d . We see that the Jacobian for this virtual system is simply A(x(t), x * ) andA(x(t), x * ) + A(x(t), x * ) ⊤ = -2P s (x(t), x * ) ⊤ P s (x(t), x * ) -2αI d ⪯ -2αI d .In other words, in view of (3), and since M (x) = I d , any two solution trajectories y 1 (t) and y 2 (t) of the virtual system satisfy∥y 1 (t) -y 2 (t)∥ 2 ≤ e -αt ∥y 1 (0) -y 2 (0)∥ 2 .Note that we can pick one trajectory to be y 1 (t) = x * for all t and we can pick y 2 (t) = x(t). Since x 0 was arbitrary, this argument establishes the claim.Clearly, the bound (12) implies exponential convergence of trajectories of the dynamical system (10) to x * . Moreover, this bound exactly characterizes equilibrium contractivity, as was defined in [10] . Notably, using the language of logarithmic norms, Theorem 33 in [10] establishes a similar result to Theorem 3 without invoking a virtual system.Note that although Theorem 3 establishes global exponential stability and equilibrium contraction, it does not establish global contractivity. Indeed, the contractivity condition ( 2) is not guaranteed to hold with M (x) = I d without further assumptions on the structure of A(x, x * ) in (10) . Remarkably, however, due to a converse contraction theorem of Giesl, Hafstein, and Mehrabinezhad, it turns out that one can construct a state-dependent M such that the dynamics are contracting. Specifically, since trajectories of the dynamical system satisfy the bound (12) , for any matrix-valued mappingEQUATIONadmits a unique solution for each x [16] . In other words, the unique M solving this matrix PDE serves as the contraction metric and will satisfy (2) with suitable choice for c. The following proposition provides the explicit solution for M in terms of the solution to the matrix PDE.Proposition 4 (Theorem 2.8 in [16] ). Let C : R d → R d×d be smooth and have C(x) = C(x) ⊤ be a positive definite matrix at each x. Then M : R d → R d×d given by the formulaEQUATIONis a contraction metric for the dynamical system (10) on any compact subset containing x * , where τ → ϕ(τ, x) is the solution to the ODE with ϕ(0, x) = x and τ → ψ(τ, x) is the matrix-valued solution toEQUATION)While it is challenging, in general, to compute the metric (15) , numerical considerations for approximating it arbitrarily closely are presented in [17] . In practice, one would select C(x) = I d and evaluate all integrals numerically. For ELCD, unless one directly needs to know the contraction metric for application purposes, it is not required to compute the contraction metric at any point during either training or inference.We remark that our parametrization for A in ( 11) is similar to the parametrization for the Jacobian of f that was presented in [6] . There are however a few key differences. Notably, A(x, x * ) may be asymmetric while the Jacobian in equation ( 3) in [6] is always symmetric. Notably, since the Jacobian in [6] is symmetric and negative definite, the vector field ẋ = f (x) is a negative gradient flow, ẋ = -∇V (x), for some strongly convex function V . On the contrary, our dynamics (10) can exhibit richer behaviors than negative gradient flows in view of the asymmetry in A(x, x * ). Additionally, since the dynamics in [6] with their parametrization can be represented as ẋ = -∇V (x) for some strongly convex V , it is guaranteed to be contracting in the identity metric, M (x) = I d [42] . On the other hand, our dynamics (10) is not necessarily contracting in the identity metric and instead is contracting in a more complex metric given in (15).Realistic dynamical systems and their flows including the handwritten trajectories found in the LASA dataset, are often highly-nonlinear and may not be represented in the form (10) or obey the bound (12) . One solution to these challenges is to transform the system to a latent, possible lower-dimensional, space and learn an ELCD in the latent space.Latent space learning is possible because contraction is invariant under differential coordinate changes. From Theorem 2 of [29] , given a dynamical system ẋ = f (x), f : R d → R d , with f satisfying (1), the system will also be contracting under the coordinate change z= ϕ(x) if ϕ : R d → R d is a smooth diffeomorphism. Specifically, if ẋ = f (x)is contracting with metric M , then the system that evolves z is contracting as well with metric given by M (z) = Dϕ(z) -⊤ M (z)Dϕ(z) -1 , where z = ϕ(x) and Dϕ is the Jacobian of the coordinate transform. In other words, We can learn vector fields that are contracting in an arbitrary metric by training a vector field f which is parametrized as (10) and using a coordinate transform ϕ.NCDS [6] , treats the coordinate transform as a Variational Autoencoder (VAE) [23] . and a model with a transform (right) when trained on data that is generated from a vector field that is contracting in a more general metric.To demonstrate the need for a coordinate transform, we consider the task of learning a vector field that fits trajectories generated by the linear system ẋ = Ax with A = -1 4 0 -1 . Clearly, this linear system cannot be represented in the form (10) with parametrization (11) . To see this fact, we observe that A + A ⊤ is not negative definite and thus no choice of P s or P a can represent this linear system. To remedy this issue, one can take the linear coordinate transform z = ϕ(x) = P x, whereP = 1 0 0 4. Then the z-dynamics readEQUATIONand now the symmetric part of P AP -1 is negative definite and thus we can find suitable choices of P s , P a . Specifically, we can take α ∈ (0, 1/2), let P a (z) = 0 0 1/2 0 z, and let P s (z) = Qz whereQ is the matrix square root of 1 -1/2 -1/2 1 -αI 2 .Note that in this toy example, asymmetry is essential to exactly represent these dynamics in the latent space, z. If we used the parametrization in [6] , we would not be able to represent these dynamics since they have an asymmetric Jacobian. We demonstrate this routine in Figure (1) . We generate two trajectories starting at (0, 2) and (0, -2), and use that data to train two ELCDs, one without and one with a learned, linear transform. Figure 1(a) shows the vector field and trajectories corresponding to an ELCD with no transform. The trajectories are forced to the center sooner than the actual data as a consequence of the bound (12) . Learning a coordinate transform allows the ELCD to learn a system that is contracting in an arbitrary metric and exhibit overshoot. This is demonstrated by the learned system in Figure 1 (b), which perfectly matches the data.There are several popular neural network diffeomorphisms including coupling layers [12, 32] , normalizing flows [40] , M-flow [6, 8] , spline flows [14] , and radial basis functions [5] .A coupling layer ϕ : R d → R d consists of a neural network θ : R k → R N with 1 < k < d and a neural network g : R N × R → R. The transform ϕ maps input x ∈ R d to y ∈ R d with the following procedure:(i) Set y i = x i for 1 ≤ i ≤ k for some 1 < k < d. (ii) Set y i = g(x 1:k , x i ) for k ≤ i ≤ dCoupling layers are invertible by doing the above process in reverse, so long as g is invertible. A coupling layer can exhibit a wide variety of behaviors depending on the choice of g.Polynomial spline curves are another diffeomorphism that are commonly used as g in coupling layers.A polynomial spline has a restricted domain which is divided into bins. A different linear, quadratic [31] , or cubic [13] polynomial covers each bin. [14] introduce rational-quadratic splines, which constructs a spline from functions that are the quotient of two quadratic polynomials.In practice, ϕ is the composition of several of transforms. Because the coupling transform only alters some coordinates, they are often used in conjunction with random permutations. If the latent space dimension is smaller than the data space dimension [6] makes the last composite function of ϕ an 'Unpad' function, which simply removes dimensions to match the latent space. The decoder ϕ is then prepended by a 'Pad' function, which concatenates the latent variable with the necessary number of 'zeros' to match the data space.Of course, if the latent dimension is smaller than that of the data dimension, then ϕ can not be bijective. However, [6] argue that as long as ϕ is injective and has range over the dataset, then f being contracting in the latent space implies that the learned dynamics are contracting on the submanifold defined by the image of ϕ -1 . In other words, NCDS cannot be globally contracting. The consequences of this are shown in Figure (2) . While NCDS can learn to admit an equilibrium point when on the submanifold, trajectories that fall off the submanifold may diverge. ELCD, in contrast, is contracting to the correct equilibrium point during all phases of training. Our task is to simultaneously learn the contracting system f (x) and the metric in which f contracts. As previously discussed, the metric is implicitly determined by ϕ(x).[6] uses a two-step method. First, they treat ϕ as a variational autoencoder and maximize the evidence lower bound (ELBO). They then fix the encoder and train the model to evolve the state in latent space. Note that the VAE objective is to make the encoded data to match a standard-normal distribution. Notably, VAE training does not encourage the encoder to transform the data to space in which the data corresponds to trajectories of a contracting system. If the data is not contracting in the transformed space, a contracting model will not be able to fully fit the data. This explains why, in our implementation of the two step-training, we are unable to reasonably learn the data. For further comparison with NCDS, we will train the encoder and model jointly. 
experiments: 
 The code for our model and data can be found here: https://github.com/seanjaffe1/ Extended-Linearized-Contracting-Dynamics. For all datasets we compare our method against NCDS [6] , Stable Deep Dynamics (SDD) [30] and Euclideanizing Flow (EFlow) [32] . See A.2 for a detailed discussion of the methods. We report the dynamic time warping distance (DTWD) [20] (see A.1) and standard deviation between predicted and data trajectories. See A.5 for model implementation details. See also the Appendix for more details on these models. We experiment with the LASA dataset [26] , which consists of 30, two-dimensional curves. We use three demonstration trajectories of each curve. As in [6] , the first few initial points of each trajectory are omitted so that only the target state has zero velocity. we stack two and four LASA trajectories together to make datasets of 4 and 8-dimensional trajectories, respectively. All data is standardized to have a mean of zero and variance of one. We use in total 10 2D curves, 6 4D curves, and 6 8D curves. Some 2D-LASA trajectories and their respective trained ELCD trajectories and vector fields are visualized in Figure (3) .We also experiment with simulated datasets. We simulate 6 trajectories of a 2,4, and 8-link pendulum (4D, 8D, and 16D respectively) each and 4 trajectories of 8D and 16D Riemannian gradient descent dynamics on a generalization of the Rosenbrock function (see Appendix A.4). Each model is trained on all trajectories of the same dimension, and then predictions are made starting from every initial point.Table 1 presents our results. ELCD performs the best in all tasks. These results shows the benefit of the increased expressiveness allowed by the skew symmetric component of our parametrization. These benefits are more apparent in the pendulum dataset. The oscillatory behavior of the pendulum dynamics is a product of complex eigenvalues in its Jacobian. Our model's skew symmetric component is what enables it to learn the pendulum dynamics so well. A sample of demonstration and learned pendulum trajectories is shown in figure (4) .The Rosenbrock dynamics are stiff and difficult to learn. In our training of SDD, we observed lack of stability and convergence, resulting in very large DTWDs. Hence, we report NaN for the results of the SDD models on the Rosenbrock dataset. While our model performs the best, there is still room for improvement. Handling such dynamics with multiple time scales is still an open challenge.Figure 4 : Phase plots of the two link-pendulum trajectories (blue) and the trajectories produced by ELCD (red). Each column is a different trajectory. The top row is the first link and the bottom row is the second link. The x-axis is angle and the y-axis is angular velocity. 
limitations: 
 Our method assumes that the underlying dynamics of the data trajectories are contracting. However, the diffeomorphism allows for a wide class of stable systems to be represented. Our method is currently implemented for trajectories that converge to the same fixed point. But, this can be overcome by manually changing the fixed point, depending on which trajectory the input data came from. Also, right now our method is limited to dynamics in R n , but we imagine that an extension to more general manifolds should be possible.As we are motivated by applications in imitation learning and robotics, in this work, we are primarily interested in problems where the underlying systems should be robustly stable, especially away from training data. For this reason, we focus on learning dynamics which are guaranteed to be globally contracting and do not address learning other types of systems, such as those with multiple fixed points, limit cycles, or chaotic attractors. We imagine that extensions of ELCD could capture these richer dynamical behaviors by leveraging weaker notions of contraction including local contraction (i.e., contractivity in the region of attraction of a stable equilibrium), k-contraction (contraction of k-dimensional bodies) [43] , transverse contraction [28] and contraction in the Hausdorff dimension [44] . Moreover, the notion of translation invariance mentioned could be studied using semicontraction theory [11] , i.e., contraction to a subspace. 
conclusions: 
 In this paper, we introduce ELCD, the first neural network-based dynamical system with global contractivity guarantees in arbitrary metrics. The main theoretical tools are extended linearization, equilibrium contraction, and a converse contraction theorem. To allow for contraction with respect to more general metrics, we use a latent space representation with dimension of the latent space equal to the dimension of the data space and train diffeomorphisms between these spaces. We highlight key advantages of ELCD compared to NCDS as introduced in [6] , including global contraction guarantees, expressible parametrization, and efficient inference. We demonstrate the performance of ELCD on high-dimensional LASA datasets and simulated multi-link pendulum and Rosenbrock dynamics. Our method shows consistent performance across all datasets.of Homeland Security, and by IBM. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the NSF or its federal agency and industry partners. 
 
 Appendix: We evaluate trajectories by comparing them against their reference trajectory using the normalized dynamic time warping distance (DTWD) [20] . For two trajectories x = {x i } i∈{1,...,t1} , x = {x i } i∈{1,...,t2} and some distance function d(•, •), let DTWD beDTWD(x, x) = 1 t 1 t1 i=1 min xj ∈x d(x i , xj ) + 1 t 2 t2 i=1 min xj ∈x d(x i , x j ) A.Our ELCD model equipped with a coupling layer-based transformation can represent any contracting dynamical system. It is known that any smooth nonlinear system with a hyperbolically stable fixed point can be exactly linearized inside its basin of attraction via a suitable diffeomorphism [25] . Effectively, our extended linearization parameterization is tasked with learning the stable linear part, while the coupling layers aim to learn and approximate this suitable diffeomorphism. As our ELCD parametrization can universally approximate any linear model, we simply need our coupling layer to universally approximate all diffeomorphisms. Indeed, [38] proved that the coupling layers are universal approximators for diffeomorphisms.In this section, we provide additional elaboration on the multi-link pendulum dataset and the Rosenbrock dataset.Multi-link pendulum: The multi-link pendulum is a simple mechanical system which is the interconnection of n rigid links under the force of gravity. We specifically assume that there is damping on each of these links, which makes the resulting dynamical system stable and almost all trajectories converge to the downright equilibrium point. In this case, the state of the pendulum is described by the n pairs (θ i , θi ), where θ i is the angle of the i-th link and θi is its angular velocity. Thus, this dynamical system is 2n dimensional.As was done in [30] , we adapt the code from http://jakevdp.github.io/blog/2017/ 03/08/triple-pendulum-chaos/ to generate data for n links.Rosenbrock datset: The classical Rosenbrock function, see [34] fEQUATIONis a nonconvex function with global minimum at (1, 1). Despite its apparent nonconvexity, it is known that a Riemannian gradient descent dynamics, is contracting with respect to a suitable Riemannian metric, see Example 3 in [42] for details. In optimization and in machine learning, the Rosenbrock function is a benchmark for the design of various optimizers and in the study of Riemannian optimization [19, 24] .For a higher-dimensional generalization of the Rosenbrock function, we considerEQUATIONwhere λ i > 0 are parameters that affect the conditioning of the function. This generalization is still nonconvex and admits several saddle points. Note that the classic Rosenbrock corresponds to n = 2, λ 1 = 1, λ 2 = 100. The unique global minimum of this generalization is at (1, 1, . . . , 1). The corresponding contracting Riemmanian gradient descent dynamics that find this global minimum areEQUATIONwhereG(x) = Dψ(x) ⊤ Dψ(x), where ψ : R n → R n is the mapping ψ(x) = ( √ λ 1 (1 - x 1 ), √ λ 2 (x 2 -x 2 1 ), . . . , √ λ n (x n -x 2 n-1 )). Note that these dynamics were designed via a generalization of the procedure proposed in [42] . To generate data, we choose four initial points and evolved them according to (24) .For the encoder ϕ, we use the composition of two coupling layers composed of rational-quadratic spline. The splines cover the range {-10, 10} with 10 bins, and linearly extrapolated outside that range. The parameters of each spline are determined by residual networks, each containing two transform blocks with a hidden dimension of 30. A permutation and linear flow layer is placed before, in the middle, and after the two quadratic spline flow layers. For ELCD, we let the latent dimension size equal the data dimension size. P a and P s are implemented as two-layer neural networks with a hidden dimension of 16.As the authors of NCDS have not yet made their code publicly available, we implemented NCDS to the best of our abilities. We used a latent dimension of size 2 for all datasets. This is in-line with the description in [6] , and necessary given the extra cost from integrating the Jacobian in higher-dimensional latent spaces. We use the same encoder ϕ as for ELCD with an added un-padding layer as the last function which removes the last d -2 dimensions.All experiments are trained with a batch size of 100, for 100 epochs, with an Adam optimizer and learning rate of 10 -3 . All computation is done on CPUs.NCDS [6] trained the model and the encoder to output x t+1 from x t using a single Euler integration step: x t+1 = ϕ -1 (ϕ(x t ) + dt * f (ϕ(x t )). Their loss was the MSE between the predicted and true x t+1 . We empirically find better performance by training our model to directly predict the velocity vector ẋt+1 = ϕ -1 J (x t )f (ϕ(x t )), where ϕ -1 J (x t ) is the Jacobian inverse of ϕ(x t ). Calculating the gradient of ϕ -1 J (x t ) efficiently required a slight alteration to the original M-flow [8] code, which we have included in our repository. Justification: See section 5Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Justification: There is no direct link between our work and a possible negative societal impact.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 


Decision: ### Summary

This paper proposes a method for learning a contractive dynamical system from demonstrations, with global contractivity guarantees. The system is flexible and parameterized by a simple neural network architecture that relies only on mild constraints. Contractivity is achieved under an implicit metric. Additionally, a diffeomorphism enables the modeling of the dynamics, as the system is learned in the latent space, where the representations help to ensure contractivity, while the mapping to the input space handles the expressivity.

$$$$

### Strengths

**Contribution**: This paper considers an important problem and addresses some limitations of the seminal NCDS model (efficiency, expressiveness, and global contraction guarantees). It introduces an innovative combination of latent space learning and extended linearization, enabling the capture of contracting dynamics through learned representations. The provided results are promising.

**Theoretical Clarity**: The paper is well-written with a clear structure (easy to follow). The theoretical guarantees and insights were appreciated by most reviewers. The concept of extended linearization is neat as well as the use of the converse contraction theorem to prove the main contractivity result.

$$$$

### Concerns

**Experiments**: Most reviewers found the experimental section to be somewhat limited, noting that it does not sufficiently demonstrate the claimed advantages of the proposed model.

**Novelty**: Concerns were raised by Reviewer HPgs, who described the proposed methodology as "incremental," noting that the benefits cannot be adequately assessed from the provided experiments.

$$$$

### Conclusion
Considering the contributions of the paper --particularly the modeling choices and theoretical results-- along with the efforts made to address the reviewers' questions and the additional experiments/extensions provided in the rebuttal, I believe this paper is a reasonable contribution to the conference and recommend acceptance as a poster. While the experimental section may be lacking, the other aspects were well-received by most reviewers, even though the approach appears somewhat incremental to previous works.

Reviews: review: paper_summary: paper_summary: This paper proposes a novel parameterization of the extended linearisation form of a dynamical system that guarantees global contractivity. Whilst the most basic form of this parameterization only ensures contractivity in some (implicitly defined) metric, a latent space version is also proposed that enables more flexible learning of dynamics that are contractive in more general metrics. Experiments are provided that show superiority over related methods.
strength_weakness: strength_weakness: The authors very quickly gloss over the motivations for learning contractive dynamics; I had to go read the NCDS paper to understand the real-world relevance of such dynamics. Let me write down below what occurred to me when reading the intro (which mostly betray my lack of experience in potential practical applications of this paper), and the authors may or may not want to modify their introduction to preempt this potential confusion.

My initial thoughts when reading the first paragraph of the intro, in particular “it is desirable to ensure [...]”:
- For one, one might actually be learning from data that comes from a dynamical system with multiple stable fixed points, so why is it desirable to enforce a single one in the learned model?
- Relatedly, why does the stability of the _learned_ dynamical system matter? I would assume that what typically matters in robotics/control is that the actual physical system implemented (and in particular, the closed control loop) be stable and robust? What does that have to do with the learned model? Is model stability important because an unstable model (where the ground truth is stable) would extrapolate very badly outside the region where there is data?

(I was thinking of a situation where deep learning is used to learn an internal model of a plant from observational data, e.g. for subsequent control purposes; in which case (i) autonomous dynamics are not necessarily very relevant and (ii) it's unclear why those dynamics would be contractive in general; now I understand that there are applications in e.g. imitation learning where it might be useful to ensure that the learned dynamics have such guarantees, but it was hard to guess from the current intro; the intro of the NCDS paper helped).
questions: questions: In the partial contraction analysis of Eq 13, it might be worth stating that x can be any function of time (as opposed to a constant). This seems critical to the eventual conclusion, where you say "we can pick y2(t) = x(t)", as this statement implicitly relies on a trajectory of the original system also being a valid trajectory of the virtual system. I'm guessing this is very standard but as someone not familiar with this type of partial contraction analysis I thought I'd mention it.

A few typos
- l.41 “enforce an exponentially decay the Lyapunov function”
- l.52 “are constrainted”
- l.53 “there are points [...]” → there may be points [...]
- l.161 “prior words” → prior works
- l.258: diffeopmorphisms
- l.287 minimize the ELBO → maximize the ELBO
- l.290 the data is corresponds
- l.456: metrc
- just a suggestion for your acronym: "ELCoD" would be easier to pronounce ("ell-cod") than "ELCD"
limitations: limitations: Could the authors please comment on whether requiring contraction behaviour (and in particular, the existence of a unique stable equilibrium) might be too constraining in some applications?

I wonder if the authors could also say a few words about the controlled case -- i.e. driven systems as opposed to autonomous systems. E.g. if the autonomous system $\dot{x} = f(x)$ obeys the guarantees provided by this work, what guarantees might be automatically inherited by the corresponding control affine system $\dot{x} = f(x) + B u$?
paper_summary: paper_summary: The paper presents a method for learning contracting representations from a dynamical system. The novelty of the method lies in learning a linear map and a coordinate transform map (Diffeomorphism), which extends the types of contracting systems that can be learned. The experiments show that the proposed methods performs better empirically than Neural Contractive Dynamical System and other approaches for learning contracting dynamics.
strength_weakness: strength_weakness: The main weakness of the paper is the experimental part. 
There are details missing and it feels "rushed", e.g., table caption unclear (see below).
Additionally, he set of tasks are very limited and don't go beyond the "usual suspects" in stable dynamical system learning (e.g. pendulum chains).

Besides the experiments, a point of limitation would be that the approach itself is by no means "groundbreaking" but rather incremental in the context of learning contracting dynamics (whether empirically it is truly groundbreaking cannot be assessed from the provided experiments)
questions: questions: The caption in Table 1 misses which metric is represented here and also what the +- means (this should not be implicitly assumed)
limitations: limitations: The paper discusses some extensions to having multiple fixed points.
However, the paper lacks a true discussion on limitations and scenarios where the approach may fail (not necessarily due to the approach but also including due to learning it)
paper_summary: paper_summary: The paper presents Extended Linearized Contracting Dynamics (ELCD), a dynamical system with neural network components that has global contraction guarantees. They demonstrate improved efficiency and performance on the trajectory fitting LASA dataset (upto 8 dim), pendulum dataset (upto 16 dim), Rosenbrock dataset (upto 16 dim).
strength_weakness: strength_weakness: 1. If the authors would like to present more of a empirical contribution, it would help to scale the experiments to something of the order of 100s of dimensions such as in common RL environments like mujoco-humanoid.
questions: questions: Please see above.
limitations: limitations: NA
paper_summary: paper_summary: The paper proposes a novel neural network architecture, Extended Linearized Contracting Dynamics (ELCD), for learning dynamical systems with global contractivity guarantees. The authors leverage the concept of extended linearization to factorize the vector field and enforce negative definiteness of its symmetric part, ensuring global exponential stability and equilibrium contraction. To extend contractivity to more general metrics, they introduce a latent space representation and learn diffeomorphisms between the data and latent spaces. The method is evaluated on various datasets, demonstrating superior performance compared to existing models, including Neural Contractive Dynamical Systems (NCDS).
strength_weakness: strength_weakness: The assumption of contractivity in the underlying dynamics might limit the applicability of the model to real-world systems where this assumption may not hold.

The current implementation is restricted to trajectories converging to a fixed point and dynamics in Euclidean space, which could be extended to more general settings. This limitation is particularly evident in the LASA handwriting dataset, where the model struggles to capture the translation invariance inherent in human handwriting. In reality, the end-point location of handwriting is not fixed, and the entire letter can be translated without affecting its validity.  Also, handwriting can often exhibit overlapping trajectories, which cannot be modeled by the current implementation. Moreover, the learned model can only generate one handwriting.

It is unclear why the joint training of the encoder and the model is possible for ELCD but not in NCDS. The authors could address this point to clarify the advantages of their approach.

The paper lacks a clear discussion of the specific types of problems where ELCD would be most applicable and where it might not be suitable. A more detailed analysis of the model's strengths and weaknesses in different scenarios would be beneficial.
questions: questions: 1. Is the ELCD model expressive enough to represent all contracting dynamical systems? If not, what are the limitations, and can they be formally characterized?

2. The paper states that "contraction is invariant under differential coordinate changes..." (L227) but also mentions that "If the data is not contracting in the transformed space..." (L291). These statements seem contradictory. Can the authors clarify this discrepancy and explain how the choice of diffeomorphism affects the contractivity of the learned dynamics?

3. The computational complexity of approximating the contraction metric is mentioned as a potential limitation. Can the authors elaborate on the scalability of the method to high-dimensional systems and suggest possible strategies to mitigate this issue?

4. The paper focuses on trajectories converging to a fixed point. How can the method be extended to handle more complex dynamical behaviors, such as limit cycles or chaotic attractors?

5. Why is joint training of the encoder and the model possible for ELCD but not in NCDS?
limitations: limitations: The limitation sections address some of the limitations of the current approach, but is insufficient. 

The current implementation of ELCD is limited to trajectories converging to a fixed point and dynamics in Euclidean space. This restricts its applicability to problems where these assumptions hold. For example, while ELCD is applied to model handwriting in the LASA dataset, it struggles to capture the translation invariance inherent in human handwriting, where the end-point location is not fixed. It would be beneficial to  discuss the specific types of problems where ELCD would be most applicable and where it might not be suitable.

Comments: comment: I thank the authors for the clarification in their response.  I'm keeping the original score.
Thank you. I will keep my score.