Learning Lyapunov-Stable Polynomial Dynamical Systems Through Imitation
Authors: Amin Abyaneh
Hsiu-Chin Lin
Venue: CoRL 2023 Poster (2023)
Keywords: Imitation learning
Safe learning
Motion planning
Dynamical system
Semidefinite programming
Robotic manipulation

Abstract: Imitation learning is a paradigm to address complex motion planning problems by learning a policy to imitate an expert's behavior. However, relying solely on the expert's data might lead to unsafe actions when the robot deviates from the demonstrated trajectories. Stability guarantees have previously been provided utilizing nonlinear dynamical systems, acting as high-level motion planners, in conjunction with the Lyapunov stability theorem. Yet, these methods are prone to inaccurate policies, high computational cost, sample inefficiency, or quasi stability when replicating complex and highly nonlinear trajectories. To mitigate this problem, we present an approach for learning a globally stable nonlinear dynamical system as a motion planning policy. We model the nonlinear dynamical system as a parametric polynomial and learn the polynomial's coefficients jointly with a Lyapunov candidate. To showcase its success, we compare our method against the state of the art in simulation and conduct real-world experiments with the Kinova Gen3 Lite manipulator arm. Our experiments demonstrate the sample efficiency and reproduction accuracy of our method for various expert trajectories, while remaining stable in the face of perturbations.

Learning Lyapunov-Stable Polynomial Dynamical Systems Through Imitation
Amin Abyanch
Hsiu-Chin Lin
Department of Electrical and Computer Engineering
School of Computer Science
McGill University
McGill University hsiu-chin, lin@cs wwgill.lll.... ca
Abstract: Imitation learning is paradigm to address complex motion planning problems by learning policy to imitate an expert's 's behavior. However, relying solely on the expert's data might lead to unsafe actions when the robot deviates from the demonstrated trajectories. Stability guarantees have previously been provided utilizing nonlinear dynamical systems, acting as high-level motion planners, conjunction with the Lyapunov stability theorem. Yet, these methods are prone inaccurate policies, high computational cost, sample inefficiency, or quasi stability when replicating complex and highly nonlinear trajectories. To mitigate this problem, we present an approach for learning globally stable nonlinear dynamical system as motion planning policy. We model the nonlinear dynamical system as a parametric polynomial and learn the polynomial's 's coefficients jointly with a Lyapunov candidate. To showease its success, we compare our method against the state of the art in simulation and conduct real-world experiments with the Kinova Gen3 Lite manipulator arm. Our experiments demonstrate the sample efficiency and reproduction accuracy of our method for various expert trajectories, while remaining stable in the face of perturbations.
Keywords: Imitation learning, Safe learning, Motion planning, Dynamical system, Semidefinite programming, Robotic manipulation
Introduction
Motion planning for robotic systems is generally regarded as decomposition desired motion into a series configurations that potentially satisfy set constraints [1]. Imitation learning tackles motion planning by imitating an expert's 's behavior to learn a planning policy [2]. To this day, only handful of imitation learning methods provide mathematical stability guarantees for their resultant policy. Stability is critical factor when deploying imitation policies in environments exposed to external perturbations. Therefore, unpredictable environments require policy that reasonably responds in unexplored regions of state space, away from original demonstrations.
Researchers have turned to autonomous dynamical systems (DS) as means learn stable motion planning policies [3, 5]. Essentially, a parametric time-invariant DS is optimized provide action (velocity) given the current state (position), while adhering to constraints that attain global Lyapunov stability. This approach leads to safety and predictability of planned trajectories, even in areas of state space without expert demonstrations. However, previous work is mostly confined to basic Lyapunov functions that adversely impact the reproduction accuracy, and require sufficiently large set of demonstrations. Others have proposed approaches focused on diffeomorphism and Riemannian geometry [6. 7, 8] and contraction theory [9], that are prone to quasi- stabiiit increased computational time, or restricted hypothesis class.
We propose a method to simultaneously learn polynomial dynamical system (PLYDS) and poly- nomial Lyapunov candidate to generate globally stable imitation policies. Polynomials, depending on
7th Conference Robot Learning 2023), Atlanta, USA.

---

Expert Demoonttatiio Dataset
Policy Formulation
Manipulater Arm
pobreciial with aquaaale croplenly
Optinization
iepiccceed
Torget
X
caribelute progoreiig
Testock
Juint Controller
kyspunay Candisaale
Actio
Pelemanii
velecty
Figure 1: Overview of the stable policy learning framework. Policy learning (left) optimizes stable polynomial DS from expert's demonstration data. This policy is then deployed (right) to plan globally stable and predictable trajectories in the entire state space.
the degree, possess an expressive power to approximate highly nonlinear systems, and polynomial regression can empirically compete with neural networks on challenging tasks [10, 11]. Unlike most neural policies, global stability can be naturally expressed with polynomials. Polynomials also enable us utilize efficient semi-definite programming [12, 13] and uum-of-qquares (SOS) optimization techniques [14, 15], and offer adaptability expert's 's demonstrations.
Our main contribution is twofold. We propose a polynomial representation of the planning policy and Lyapunov candidate function, coupled with concrete mathematical stability certification for precise and safe replication of the expert's demonstrations, as depicted in Figure define a regularized semi-definite optimization problem to jointly learn the DS and the Lyapunov candidate with higher flexibility and precision. We compare the reproduction accuracy of PLYDS with alternatives in the literature and evaluate the performance in both simulation and real robotic systems.
Then, we
2 Background and Notation
Consider system operating in a state-space X' C RR, e.g., robot in its task- or configuration-space. The system can execute actions in A C R", for instance, velocity or torque commands, leading to state evolution. We denote the state variable with x - [x] I2 xn)² E X, and consider the action variable to be the state's derivative E A. Within this space, our goal is to learn an imitation policy through a dataset of experts' state-action pairs, referred as trajectories.
Let Na N be the number of trajectories demonstrated by the expert. Each trajectory contains N, E N state-action pairs. The dataset of expert trajectories stacks all state-action pairs, defined as:
D a (x)),, (s)) d (1...., Na), E (1, (.....,N.)..
where (s), (s)) the dataset entry corresponding the s-th sample of the d-th demonstrated trajectory. The dataset D holds N Na samples. We assume that the trajectories contain the same sample size (N.), share common target (x" E X), and have zero velocity at the target, i.c., x²(N.) x" and s Definition 2.1. (Dynamical Systems). The mapping between the state and the action in each sample can be modelled with time-imeen invariant autonomous dynamical system (DS), denoted by:
0 for all trajectories d E (1....,Nd)
x =
f(x), x
In Equation (2), f is an ordinary differential equation for the true underlying DS. The term E IR" captures measurement and recording noise of expert's demonstrations. We assume that e embedded in the estimated DS, f. and eliminate the need for modeling its distribution. Following [3], we aim at learning a noise-free estimation of (x), denoted by f(x). One can view f(x) in Equation (2) as policy that maps states to actions for reproducing the demonstrated trajectories in the state-space. For instance, when the robot is located in xa E X, the policy yields action Xn = f(xo), which can be passed to the robot's 's velocity controller.

---

The estimated DS in Equation (2), fxx, is globally asymptotically stable (GAS) around an equilib- rium point x if and only if for every initial state, x + x" as the system evolves and time goes to infinity [16]. A popular tool to study the GAS property of a DS is the Lyapunov stability theorem. According to this theorem, a DS exhibits GAS there exists a positive-defintte function v w: X - R, known as Lyapunov potential function (LPF), such that i(x) < 0 for all x * x" and i(xx) 0. To ensure GAS for fxx), we simultaneously learn the policy, f. and the LPF, v.
3 Related work
Extensive research is conducted on imitation learning and its applications in robotic motion planning for variety of tasks. Existing efforts can be divided into the following predominant research tracks.
Dynamical systems for motion planning. Dynamical systems have proved to effectively counter autonomous motion planning problems by proposing a time-invariant policy [17]. Traditional methods of encoding trajectories are based on spline decomposition [18], Gaussian process regression [19], or unstable dynamical systems [20, 21]. They either lack robustness because of time-variance or fail to provide GAS, SEDS [3] is the first attempt learn stable planning policies. However, its performance declines when applied to highly nonlinear expert trajectories. Most notably, it suffers from trajectories where the distance to the target is not monotonically decreasing. The intrinsic limitation of SEDS comes from the choice of simple Lyapunov function. Follow-up research introduces more complex Lyapunov candidates to stably mimic nonlinear trajectories [4, 22], but are still restricted in representing the Lyapunov candidate. Others have tried to tackle SEDS limitations through diffeomorphic transformations and Riemannian geometry [6, 8, 7] that yield quasi-stable planners for some trajectories, and contraction theory [9] that restricts the class metrics to make the optimization tractable. Lastly, most improvements the original SEDS still the Gaussian mixture model formulation, that is vulnerable in presence of limited expert demonstrations.
Imitation learning. Recent imitation learning developments can be applied to motion planning tasks with minimal modifications, since motion planning can be achieved by finding a (not necessarily stable) policy in the robot's 's spske from the expert's behavior. For instance, GAII [23] introduces an adversarial imitation learning approach that directly optimizes the expert's policy, but requires large set of expert's 's data (low sample efficiency) and extensive training iterations. The growing interest in neural policies has also led to the development of end-to-end autonomous driving [24] and behavioral cloning [25, 26, 27] methods. Nevertheless, they generally lack GAS, and it is unclear whether the robot can recover from perturbations. The same drawbacks exist with apprenticeship learning approaches, such as Abbeel and Ng [28] and inverse reinforcement learning, such as Ziebart et al. [29], and the computational demand is even higher for the latter.
Stability in neural dynamical systems. Methods such as [30, 31] represent the dynamics with Neural Network, and propose the joint training of dynamics and Lyapunov function to guarantee the stability. Though theoretically sound, these methods have only been applied to rather simple settings and require large demonstration sets. Neural Lyapunov methods [32, 33, 34] promise a data driven and potentially stable approach to control and model nonlinear dynamies, but lack global stability. Methods such as [35] are also not stable-by-design and the dynamical system lacks autonomy.
4 Methodology
We instantiate the policy and the corresponding LPF candidate, f and v, with polynomials in Section 4.1 and Section 4.2, respectively. This allows us to accurately imitate expert's behavior, while providing formal GAS guarantees. Subsequently, we formulate tractable optimization problem for jointly learning the policy and the LPF in Section 4.3.
4.1 Dynamical system policy formulation
We need to approximate the unknown underlying DS in Equation (2) to discover the mapping between states and actions from expert's behavior. To this end, we opt to model the policy with
3

---

parametric polynomial. The representative power of polynomials was originally established through the Weierstrass approximation theorem, stating that every continuous function defined on a closed interval can be approximated with desired precision by a polynomial. This idea is fortified by recent studies, such as [10, 11], that compare polynomials to neural networks on a variety of tasks.
Definition 4.1. (Polynomial Dynamical Systems). Polynomial Dynamical System (PLYDS) is polynomial approximation of the policy in Equation (2), and is expressed as,
x = f(x; P) - 1b Pibxa by P2bx, ... "x.aaa
(3)
where bx,a [1 (x (x²) (x is the element-wise k-th power of xT Every row i of f is a polynomial of degree 2o, fi(x; P:) = by P.bx,a where P; E gon+1 and Sh a {S RA*KIST =S}. The matrix P P gan"+" encapsulates the block-diagonal form of all P; matrices.
(aaf) is the polynomial basis of degree a E N, and
Below, we present an example to show how PLYDS, as defined in Definition 4.1, captures nonlinear time-invariant policies. One can further complicate the policy by increasing o, which in turn produces a larger basis vector and a more flexible polynomial.
Example 4.1.1. A second-order polynomial representation of a one-dimensional DS is:
X = f(x; P) [bT P1 bx.a] =
[1
x]
Poo
P01
1
=oo== + (poi + P1O) + poo,
Pun
Pn1
where a = 1, bx,a = [1 x."". Note how P can be symmetric without loss of generality.
4.2 Global stability guarantees for polynomial dynamical systems
As explained in Section 4. a polynomial policy allows for accurately imitating the expert's demon- strations. Yet, there is no formal GAS guarantee that the robot will ultimately converge the target in the face of perturbations, deflecting it from the expert's trajectories. Owing to the Lyapunov stability theorem, finding an LPF that meets the criteria in Section 2 ensures the desired stability [36].
The major challenge lies in learning an LPF, e, which is a positive definite function with negative gradient. We tackle this by confining to the class of polynomial LPF candidates. Definition 4.2. (Polynomial Lyapunov Candidate). A multidimensional polynomial LPF is given by,
v(x: Q) A 2 Q2bxa b v: X - + R
(4)
where B E N is the polynomial basis degree. Each row is defined by ve(x; Q₂) = bx Q bx, u
X' + R, and can be viewed as a scalar Lyapunov function. The parameters malrix, Q E gBn +n is a dlako-diagonal of all Q2 E gSn+1 matrices.
Definition 4.2 introduces non-conventional LPF candidate. Rather than considering a single LPF, we designate a distinct polynomial LPF for each dimension of the state space and stack them into vxx; Q). This characterization, known vector Lyapunov function [37], is less restrictive and enables the policy and LPF to be learned moreindependently for each dimension of the state space.
We highlight that the GAS of the policy in each dimension, fx(;; P,), implies the GAS of the entire policy, f(x; P). Proposition 4.3 establishes a link between the policy stability in each row and the global stability of the multidimensional policy.
Proposition 4.3, Assuming each pair (f;(x; Pt), (x(:;; Q:)) individually satisfies the GAS condi- tions. Then, the sum ê = Σ vx(x; Q:) yields valid standard Lyapunov function for f(x; P), proving that the policy satisfies GAS conditions. The proof is given in Appendix A.I.
The formulation of the policy and the LPF as multidimensional polynomials empowers us to leverage tools from sum-of-squares (SOS) [15, 38]. The SOS approach boils the Lyapunov GAS conditions down to verifying -esitive-eefnniteness of a set of specified matrices. The next two lemmas illustrate the SOS formulation of Lyapunov stability conditions.

---

Lemma 4.4. The first Lyapunov stability criterion, vx(x; Q.) = O. is satisfied for each i e (1..... n) Q / O and Q E 30n+1. The proof is outlined in Appendix A.2.
Lemma 4.5. The second Lyapunov criterion, at (x; Q₂) < O, is fulfiled for each i (I....,n]
if theeee exists a symmetric matrix G; - 0 and G; E gla++)))++ such that:
① (x; Q.) @v;(x; Q:) θx 0x;(x; Qi) f(x: - bx,4+ b Gebx, 01
(5)
Ox al
8x
where a + 3 is the basis degree. The matrix G; is acquired by polynomial coefficient matching, and depends on and Q We summarize this dependence with the function G(P, Q) = G, where G symbolizes the block diagonll form of all G, matrices. The proof outlined in Appendix A.3.
Finally, with the necessary tools our disposal, stability of the policy and finding SOS polynomials in Theorem 4.6. This theorem serves as the fundamental basis for the subsequent policy optimization process.
establish the connection between the global
Theorem 4.6. A polynomial DS policy, f(x; P), is GAS if the following conditions are satisfied:
(c) G(P, Q) = G. The proof is straightforward and is sketched in Appendix A.4.
(a Q = 0, (b) G <0,, O,
4.3 Joint optimization problem
At this stage, we have established polynomial representations for both the policy and the LPF, along with a firm connection that confirms global stability. Now, we develop an objective function using the Mean-Se -quan-Suarred Error (MSE) cost with the Elastic Net Regularization [39]. The MSE is calculated between the policy output and the expert's actions across demonstrated trajectories, and it solely depends on the policy parameters. Essentially, this problem entails regularized polynomial regression to minimize the imitation MSE to expert's 's demonstrations, subject to the existence of an LPF that satisfies the Lyapunov conditions. The optimization problem is framed as:
1 Na
min J(P) Q.G.P
P) - -x)
A1 P + Ao P
2N,
s.1. (a) Q G+) (P) = G (d) Q = Q", G= GT P-
where |l111 and ||- denote the first and Frobenius norms, and \ X2 E R represent the regulariza- tion coefficients. Equation 7 is a semi-definite programming with nonlinear cost function [40, 38], and can be solved using standard semi-definite solvers [41, 42]. Semi- cemi-definite programming facilitates optimization over the convex cone symmetric, positive semi- cemi-defintt matrices or affine subsets. Note that (c) can cause the optimization to become non-convex. To alleviate this, we employ SDP relaxations [12], iterative methods based on initial guess of the Q matrix, and ultimately sequential quadratic programming (SQP) [43].
Notice that the negative-definite constraints can be transformed to a semi-definite constraints by a constant shift. Furthermore, the Lyapunov conditions restrict the gradient to nonzero values away from the origin. The restriction ensures that the LPF has only one global minimum at the target.
5 Experiments
We employ two motion planning datasets for experiments. Our primary data comes from the widely recognized LASA Handwriting Motion Dataset [44], which comprises data recorded from handwritten trajectories. The second dataset contains expert demonstrations collected through teleoperating a robotic arm on realistic manipulation tasks. Details about both datasets can be found in Appendix I.
5.1 Evaluation
For evaluation purposes, we apply PLYDS and the baselines to the dataset (Figure 2a) and evaluate the performance of policy rollouts in PyBullet simulation (Figure 2b) before deploying safe policies
5

---

bxpert Ocmorstration Paicy Watout
KINDYA
9
(a)
(b)
(c)
Figure 2: Overview of the evaluation sequence: (a) learning from demonstrated data, (b) numerical evaluation in simulation and (c) deployed in real-world Gen3 Lite manipulator.
PDDDS SEDS
LPV GAI
400
0.5
8000
Computationn
200
0.
Angle
Motion Ang
Mobion Shapes
(a)
(b)
Figure 3: Comparison of (a) mean and standard deviation of reproduction MSE and (b) computation time to designated imitation learning methods. PLYDS performs reasonably well in terms of accuracy and is even more promising in terms computational cost.
onto manipulator (Figure 2c). In all experiments, randomly split the demonstrated trajectories in the dataset into train and test sets. The policy learning stage, introduced in Equation (7), is carried out on the training data. The learned policy is subsequently evaluated by calculating the MSE between Ntaet
the policy predictions and the ground truth in the test data, 2N N,
Ee Σ (f('((s)) P)
x (s)) Recall that the policy output, x, is the velocity passed to the robot's low- controller. We repeat this procedure over 20 different random seeds, and report the average and standard deviation.
We compare the accuracy of our approach existing baselines. Primarily, compare against Stable Estimator of Dynamical Systems (SEDS) [3], Linear Parameter Varying Dynamical Systems (LPV-DS) [22]. and Stable Dynamical System learning using Euclidcanizing Flows (SDS-EF) as methods that ensure GAS. We also compare our method to Behavioral Cloning (BC) [26], and Generative Adversarial Imitation Learning (GAIL) [23] to highlight the importance global stability. Note that among these, BC and GAIL. do not provide mathematical stability guarantees, but the results could provide further comparison ground for the accuracy and computation time, The implementation details, hyperparameters and architecture are discussed in Appendix B.2 and Appendix B.3.
5.2 Handwriting dataset
We compare the learned policies of PLYDS to the baselines on eight designated motions. The outcome of these experiments is reported in Figure 3. Despite stability guarantees, the overall accuracy is better among stable imitation learning methods compared to unstable neural approaches.
To analyze GAS, we visualize the learned policies of all methods as streamlines. Figure 4 illustrates the policy rollouts for N-Shaped motion of the handwriting dataset. Each sub-figure represents trained policy illustrated with gray streamlines. It is evident that SEDS and PLYDS maintain GAS,

---

Polbyy
Ieprodicet
larget
102
PLYOS
SEDSS
S08-E3
-00
Figure 4: Policy rollout for N-Shaped motion of the handwriting dataset. Each figure represents trained policy (gray) and rollouts (red) learned from demonstrations (blue). Note the stability issues with GAIL and SDS-EF, where some streamlines fail to reliably converge to the target.
Nesrodazod
Eipet Demsmntrat ore
-03
PLYDS
SEDS
SDS-EF
Figure 5: Policy rollout for Sine-Shaped motion (blue) the handwriting dataset, with access to only one expert demonstration. Each figure represents a trained policy (gray) and one rollout (red) learned from one demonstration (blue). Methods requiring large datasets for clustering, such as SEDS and LPV-DS, exhibit inaccurate and unsteady performance.
while GAIL and SDS-EF fail demonstrate converging trajectories for the entire state space. The same pattern persists for other motions depicted in Appendix C.I.
Finally, we examine the sample efficiency of our method by reducing the input data to one demon- strated trajectory. From Figure 5, we can see that PLYDS learns stable policy with such limited training samples, while the baselines generate trajectories which diverges from expert data.
So far in this section, the policy and the LPF polynomial degrees were set to a and 3 2. To understand the way in which the complexity of polynomials affects the overall performance, repeated the same experiments with degrees of = 4, 6, and 8, and present the result in Appendix C.2. We observe that higher complexity leads to improved precision, if not halted by overfitting or stability sacrifice. Moreover, we study different LPF complexities in Appendix C.3, evaluate the performance PLYDS with noisy demonstrations Appendix C.4, and further investigate the computational times in Appendix C.5.
5.3 Manipulation tasks
To conduct real-world trials, we collect second of expert demonstrations through teleoperating Kinova Gen3 Lite, manipulator with six degrees of freedom. This new dataset holds three distinct motions: (a) root-parabola, (b) standard pick and place, and (c) prolonged-sine, which represent exemplary nonlinear trajectories (Figure 6). Additional details are available in Appendix B.
The performance of all methods is summarized in Table 1, where PLYDS often outperforms other baselines. Next, the learned policy of PLYDS transferred the physical arm (Figure 6) and successfully imitates the introduced manipulation tasks. We also start the robot regions that further away from the demonstrations and introduce perturbations by randomly pushing the robot to reveal the inherent GAS of PLYDS. As expected, PLYDS manage to successfully recover to the goal.

---

it
(a)
(b)
(c)
Figure 6: Manipulation tasks: (a) root- poot-parabola, (b) standard pick and place, and (c) prolonged
Expert Motion
Prolonged Sine
Root Parabola
Fick-and-Place
Computational Time
SEDS [3]
0.234 + 0.015
0.132 + 0.023
0.094 + 0.012
277.02 + 13.60
BC [26]
1.650 + 0.133
0.931 + 0.078
0.725 + 0.133
38.93 = 9.11
GAIL [23]
2.322 + 0.098
1.322 +: 0.094
0.663 = 0.098
143.15 = 8.68
SDS-EF [8]
0.152 0.023
0.094 + 0.012
715.62 + 18.79
LPV-DS + P-QLF [22]
0.234 0.015
0.152 0.023
0.094
+ 0.012
334.55 + 25.74
PLYDS (ours)
0.111 + 0.007
0.176 + 0.015
0.021 = 0.003
21.37 + 1.52
Table 1: Policy rollout reproduction MSE and computational time in PyBullet.
6 Conclusion and Limitations
We introduced an approach that aims to learn globally stable nonlinear policies represented by polynomial dynamical systems. We employ the learned policies for motion planning based on imitating expert demonstrations. Our approach jointly learns polynomial policy along with parametric Lyapunov candidate that verifies global asymptotic stability by design. The resulting DS is utilized as a motion planning policy, guiding robots to stably imitate the expert's behavior. A comprehensive experimental evaluation is presented in real-worid and simulation, where the method is compared against prominent imitation learning baselines.
I Limitations. A limitation of SOS that the set of non-negative polynomials is larger than the ones ex- pressed as SOS [45]. Though rare in motion planning tasks, this implies that finding Lyapunov candi- date could be difficult, especially with simultaneous search for a suitable dynamical system. Lasserre hierarchy and SOS extensions [46] can search in broader class of LPF candidates and tackle this is- sue. Another limitation occurs when finding highly complex policies that lead to a violation of stability guarantees. This often happens when the regularization coefficients or the optimization tolerance are not set properly. We discuss this trade-off between stability and accuracy in Appendix C.2. Further, the computation complexity of PLYDS is feasible with reasonable choice of polynomial degrees. Higher degrees are computationally demanding, but are often unnecessary in normal motion planning tasks.
Future work. Future work includes incorporating more elaborate safety criteria, such as control bar- rier functions [47] or real-time obstacle avoidance, into our learning objectives. Plus, applications of our method in SE(3) planning, or other higheerdimensional spaces, such as configuration space of ma- nipulator robots, may be further investigated. Vector Lyapunov functions and adaptable complexity of polynomials can pave the way for such applications, as they assuage major computational challenges.
Video, Codebase, and Reproducibility
The codebase, video supplements, etc. related to this project are available on our Git repository ! Reproducing the experiments is as straightforward installing dependent software packages, and running a Unix commands in README files.
8

---

Acknowledgments
This work sponsored by NSERC Discovery Grant. We also appreciate the thoughtful reviewers' comments which helped us enhance the paper, particularly the experiments.
References
[1] J.-C. Latombe. Robot motion planning, volume 124. Springer Science & Business Media, 2012.
[2] A. Hussein, M. M. Gaber, E. Elyan, and C. Jayne. Imitation learning: A survey of learning methods. ACM Computing Surveys (CSUR), 50(2):1-35, 2017.
[3] S. M. Khansari-Zadch and A. Billard. Learning stable nonlinear dynamical systems with gaussian mixture models. IEEE Transactions Robotics, 27(5):943 957, 2011.
[4] M. Khansari-Zadeh and A. Billard, Learning control Lyapunov function to ensure stability of dynamical system-based robot reaching motions. Robotics and Autonomous Systems, 62(6): 752-765, 2014.
[5] N. Figueroa and A. Billard. physically-consistent bayesian non-parametric mixture model for dynamical system learning. In A. Billard, A. Dragan, J. Peters, and J. Morimoto, editors, Proceedings of The 2nd Conference Robot Learning, volume 87 of Proceedings of Machine Learning Research, pages 927-946. PMLR, Oct 2018. URL https::/proceeddigg. mlr. press//87//fiuurool188 html..
[6] K. Neumann and J. J. Steil. Learning robot motions with stable dynamical systems under diffeomorphic transformations, Robotics and Autonomous Systems, 70:1-15, 2015.
[7] D. Ratliff, J. Issac, D. Kappler, S. Birchfield, and D. Fox. Riemannian motion policies. arXiv preprint arXiv: 002554444. 2018.
[8] M. A. Rana, A. Li, D. Fox, B. Boots, Ramos, and N. RatlifT. Euclideanizing flows: Dif- feomorphic reduction for learning stable dynamical systems. In Learning for Dynamics and Control, pages 630-639, PMLR, 2020.
[9] H. Ravichandar, I. Salehi, and A. Dani. Learning partially contracting dynamical systems from demonstrations. In Conference on Robot Learning, pages 369-378. PMLR, 2017.
[10] X. Cheng, Khomtchouk, N. Matloff, and P. Mohanty, Polynomial regression as an alternative neural nets. arXiv preprint 1rXiv::8806.0m
2018.
[11] Morala, J. Cifuentes, R. E. Lillo, and Ucar. Towards a mathematical framework to inform neural network modelling via polynomial regression. Neural Networks, 142:57-72. 2021.
[12] L. Vandenberghe and S. Boyd. Semidefinite programming. SIAM review, 38(1):49-995 1996.
[13] H. Wolkowicz, R. Saigal, and L. Vandenberghe. Handbook of semidefinite programming: theory, algorithms, and applications, volume 27. Springer Science Business Media, 2012.
[14] M. Laurent. Sums of squares, moment matrices and optimization over polynomials. In Emerging applications of algebraic geometry, pages 157-270. Springer, 2009.
[15] P. A. Parrilo. Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization. California Institute of Technology, 2000.
[16] R. I.. Devaney. An introduction to chaotic dynamical systems. CRC press, 2021.
[17] A. Billard, S. Calinon, R. Dillmann, and S. Schaal. Survey: Robot programming by demonstra- tion. Technical report, Springrer, 2008.

---

[18] J. Aleotti and S. Caselli. Robust trajectory learning and approximation for robot programming by demonstration. Robotics and Autonomous Systems, 54(5):409-44 2006.
[19] M. Muhlig, M. Gienger, S. Hellbach, Steil, and C. Goerick. Task-level imitation learning using variance-based movement optimization. In 2009 IEEE International Conference Robotics and Automation, pages 1177-1184. IEEE, 2009.
[20] M. Hersch, F. Guenter, S. Calinon, and A. Billard. Dynamical system modulation for robot learning via kinesthetic demonstrations. IEEE Transactions Robotics, 24(6):1433-1467, 2008.
[21] Schaal. Scalable locally weighted statistical techniques for real time robot learning. Applied Intelligence-Special issue on Scalable Robotic Applications of Neural Networks, 17:49-60, 2002.
[22] N. B. Figueroa Fernandez and A. Billard. A physically-consiiten bayesian non-parametric mixture model for dynamical system learning. Technical report, EPFL, 2018.
[23] J. Ho and S. Ermon. Generative adversarial imitation learning. Advances in neural information processing systems, 29, 2016.
[24] M. Bojarski. D. Del Testa, Dworakowski, B. Firner, B. Flepp. P. Goyal, L. D. Jackel, M. Monfort, U. Muller, J. Zhang, al. End to end learning for self-driving cars. arXiv preprint arXiv:160.007316 2016.
[25] Torabi, G. Warnell, and P. Stone. Behavioral cloning from observation. Proceedings of the 27th International Joint Conference on Artificial Intelligence, pages 4950-4951 4957, 2018.
[26] D. A. Pomerleau. Alvinn: An autonomous land vehicle in neural network. In Advances in neural information processing systems, volume I, pages 305-313, 1988.
[27] P. Florence, C. Lynch, A. Zeng, O. A. Ramirez, A. Wahid, L. Downs, A. Wong, J. Lee, I. Mordatch, and J. Tompson. Implicit behavioral cloning. In Conference on Robot Learning. pages 158-168. PMLR, 2022.
(28) Abbeel and A. Y. Ng. Apprenticeship learning via inverse reinforcement learning. In Proceedings of the twenty- international conference Machine learning, page 1, 2004.
[29] B. D. Ziebart, A. L. Maas, J. A. Bagnell, A. K. Dey, et al. Maximum entropy inverse reinforce- ment learning. Aaai, volume 8, pages 1433-1438. Chicago, II, USA, 2008
[30] J. Z. Kolter and G. Manek. Learning stable deep dynamics models. Advances in neural information processing systems, 32, 2019.
[31] N. Lawrence, P. Loewen, M. Forbes, Backstrom, and B. Gopaluni. Almost surely stable deep dynamics. Advances in Neural Information Processing Systems,
2020.
[32] M. Richards, Berkenkamp, and A. Krause. The Lyapunov neural network: Adaptive stability certification for safe learning of dynamical systems. In Conference on Robot Learning. pages 466-476, PMLR, 2018.
[33] H. Dai, B. Landry, I.. Yang, M. Pavone, and R. Tedrake. Lyapunov-stable neural-network control. In Proceedings of Robotics: Science and Systems, July . doi: 1tw.:1...607//.S.... 200.:/XS...22......0.....
[34] A. Coulombe and H.-C. Lin. Generating stable and collision-fres policies through Lyapunov function learning. International Conference on Robotics and Automation, pages 3037-3043, 2023.
[35] Y.-C. Chang, N. Roohi, and S. Gao. Neural Lyapumov control. Advances in neural information processing systems, 32, 2019.
10

---

[36] A. Papachristodoulou and S. Prajna. On the construction of Lyapunov functions using the sum of squares decomposition. In Proceedings of IEEE Conference on Decision and Control, volume 3, pages 3482-3487. IEEE, 2002.
[37] V. Lakshmikantham, V. M. Matrosov, and S. Sivasundaram. Vector Lyapunov functions and stability analysis of nonlinear systems, volume 63. Springer Science & Business Media, 2013.
[38] G. Blekherman, P. A. Parrilo, and R. R. Thomas. Semidefinite optimization and convex algebraic geometry. SIAM, 2012.
[39] H. Zou and T. Hastie. Regularization and variable selection via the elastic net. Journal of the royal statistical society: series B (statistical methodology), 67(2):301-320, 2005.
[40] S. Burer and R. D. Monteiro. A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization. Mathematical Programming, 95(2):329-33) 357,-3577777 2003.
[41] M. Bazaraa, H. D. Sherali, and C. M. Shetty. Nonlinear programming: theory and algorithms. John Wiley & Sons, 2013.
[42] M. ApS. MOSEK Optimizer API for Python 10.0.46, 2021. URL http://docs http:/occc.osskk....oo mosek, com/ wotoooooo/iineeeeeeoeooooooooooo iinpex/iiixxx.eee html.
[43] P. T. Boggs and J. W. Tolle. Sequential quadratic programming. Acta numerica, 4:1-51, 1995.
[44] S. M. Khansari-Zadeh and A. Billard. Lasa Handwriting Dataset. https://cs. http///ossssaaoord....o htm1#SEDS reference, 2011.
edu/ httpl/pppppkaaaallllii
[45] G. Blekherman. There are significantly more nonegative polynomials than sums of squares. Israel Journal of Mathematics, 153:355-3800 2006.
[46] A. Papachristodoulou and S. Prajna. On the construction of Lyapunov functions using the sum of squares decomposition. In Proceedings of the s1 IEEE Conference on Decision and Control, 2002., volume 3, pages 1482-3487 vol.3, 2002. doi: 4o111..109/CCCC222..1144
[47] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada. Control barrier functions: Theory and applications. In 18th European control conference (ECC), pages 3420-3431. IEEE, 2019.
[48] B. O' O'Donoghue, E. Chu, N. Parikh, and S. Boyd. Conic optimization via operator splitting and homogeneous self-dual embedding. Journal of Optimization Theory and Applications, 169(3); 1062-1068, June 2016. http:httt///ttfoooee www.eedod//pppprrs.eeeeeeeeeee
html.
[49] Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, Reddy, D. Cournapeau, E. Burovski. Peterson, W. Weckesser, J. Bright, S. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, i. Polat, Y. Feng, E. W. Moore, J. VanderPlas. D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris. A. M. Archibald, A. H. Ribeiro. F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272, 2020. 1oi:10. 990N441199-009-088--2
[50] S. G. Nersesov, W. M. Haddad, and Q Hui. Finite-time stabilization of nonlinear dynamical systems via control vector Iyаpunov functions. Journal of the Franklin Institute, 345(7):819- 837, 2008

---

A Mathematical Proofs
Due to space limitations and to maintain coherency, proofs for propositions, lemmas, and theorems are presented this section in the order in which they appeared the main text.
A.1 Proof of Proposition 4.3
Assuming each pair (fi(x; P.), 0x(x; Q.)) individually satisfies the GAS conditions. Then, the sum
ē Ei=1 vx(x; Q.) yields a valid standard Lyapunov function f(x: P), proving that the policy satisfies GAS conditions.
Proof. As vx(x; Q:) is an LPF candidate for fi(x; Pi), both the first and second Lyapunov conditions must be satisfied, i..., vi E {I,...,m}:
(b) (x:(x; Qi)
(a) (x: Qi) ビ 0, Vx E E
x0, Vx E X.
Øt
Define the sum of elements o(x; Q) - Σラ (x Q:). We show that o(x; Q) satisfies both Lyapunov global stability conditions:
Q.)
- 0 (a) 2n (x: Q1) ... + Ux(x: Qn) 0(x; Q) = 0,
a0 x
>Σ 1-1 U/ (x: Q.)
Ov;(x; Qi)
(x;(x: Q:)
i8
Øt 0t < 0 (b). 0
Öt
Ot
A.2 Proof of Lemma 4.4
The first Lyapunov stability criterion, (x(:; Qi) ž O, is satisfied for each i E (1....,nn) n) Q ł 0 and Q E g3n+1
Proof. Considering that (x; Q,) = bT Q, and Q: is not singular, we can perform Cholesky factorization on the parameters' matrix Q. The result Q L L,, and the positivity u; Q.) comes from,
(x; Q by Q,bx, L L;bx, (L,bx.3)* (L,bx,)) || 2 0,
that represents (x; Q₂) as an SOS and therefore achieves the first Lyapunov condition.
A.3 Proof of Lemma 4.5
The second Lyapunov criterion, 07 (x; Qi) -0 is fulfilled for each i (I....,m) if there exists symmetric matrix G; - 0 and G, E gla++))+11 such that:
① (x(:; Qu) av;(x; Qi) 0x Ov;(x; Qi) f(x: P) P)= (0)=b2 bT.04. G,bx,
(8)
0t"
Ox Øt
8x
where a B the basis degree. The matrix G; is acquired by polynomial coefficient matching. and depends on P and Q; We summarize this dependence for all E (I,.... with the function G(P,C = G. where G symbolizes the block- bock-diagonall form of all G; matrices.
Proof. We know that LPF rows are denoted by vs vx(x; Q.) = b% Q.bx.a. Hence, we write the second Lyapumov condition by taking the derivative of each row:
(x;(x; Q:) ix;(x; Qe) d.
(x;(x; Q.) dr
0v,(x: Qi) Ox Ox Øt
+
Ot dr
Öt
Ox2
Ot
0z, [x: P₂)
Sv, Ov,(x; Qi)
n
du; (x; Qi) 1x(x; P,)
Ø1
Oz,
a G Q b
[b P (x,a]
dx,
12

---

Within the last summation, both the derivative of Lyapunov function and the policy are polynomials. The idea is that their multiplication could also be written as an SOS polynomial the parameters P: and Q are chosen carefully. For this polynomial, we define new basis bx,a++3 and G, E gla++3)++1 Note that the degree of this basis calculated by [(23 - 1) + (2x) + I]. which is the rounded-up degree of the above multiplication term.
Next, we match polynomial coefficients on both sides, yielding G; parameters as function of both P and G i.e..
0 x.9
Q.b..s) [b Pybxal
by xx.0+ G;bx,++3 Coef ficienntsts
Matching
Oig
+ +G (P..Q.)
We summarize the same relationship for all G; matrices, and call the resulting function G. Hence, the second condition can be represented by G G(P, Q) and G O, and be viewed SOS. D
A.4 Proof of Theorem 4.6
Assuming the polynomial representation of nonlinear autonomous dynamical system (Defini- tion 4.1),
f(x; P) = [b Pi bx, by P2bx.
,, bx,
existence a corresponding polynomial Lyapunov function (Definition 4.2), b Q2bx, bx bx
v(x; Q) = | Qibx,a
guarantees the asymptotic global stability of the policy, the following conditions are satisfied:
aa Q 0, (b) G 0,
((, Q) G.
Proof. The proof is straightforward and follows both Lemma 4.4 and Lemma 4.5. We know that each partial DS fi(x; P:) is stable if the corresponding parameterized LPF satisfies (a), (b), and (c), where G is an affine function found in Lemma 4.5 by polynomial coefficient matching. Since each [x; P²) explains the derivative along one of the orthogonal basis of f(x; P), their individual global stability equivalent to the stability of the entire system. In other words,
Va; E D()((x Pd). lim 24 x
1
lim x lim t-000
lim D2 1-000
lim zn =xx t00
Another proof can be provided using the LPF introduced in Proposition 4.3 as a Lyapunov candidate for the whole system. Both proofs equally validate the stability of the polynomial DS.
B Experiment Setup and Details
Enclosed in this section are detailed descriptions of our experiment setup, main software packages, and datasets. Due to space limitations, crucial details from the experiments explained here. Even though reading the section is not necessary to understand the paper, it provides useful insight into our setup and can aid reproducibility and future research.
B.1 Datasets
Handwriting dataset. The LASA Handwriting Dataset, partly depicted in Figure 7, is a collection of 2D handwriting motions recorded from Tablet-PC and by user's input. The dataset includes 30 human handwriting motions, where each motion represents desired pattern. For each pattern, there seven demonstrations recorded, with each demonstration starting from slightly different (but fairly close) initial position but ending at the same final point. These demonstrations may intersect
13

---

Target
Training Demoe
G
Angle
LINC MSTS
Sine
Figure 7: Plots of handwriting dataset motions used in our experiments. We select representative subset motions for baselining to keep the experiments computationally feasible. Each plot shows demonstrations with 1000 recorded samples per each. Notice that the time indexing is included in the dataset, but irrelevant our work learn time-invariant policies.
with each other. Out of the 30 motions, 26 correspond to single pattern, while the remaining four motions include multiple patterns, referred to Multi Models. In all the handwriting motions, the target position is defined as (0, O), without loss of generality. The dataset provides the following features:
Position (2 x 1000) representing the motion in 2D space. The first row corresponds to the x-axis, and the second row corresponds the y-axis Cartesian coordinates.
Time (I x 1000) being the time-stamp for each data point in the motion. We do not use this property, as our proposed method generates time-invariant policies.
Velocity (2 x 1000) representing the velocity corresponding to each position. We use this feature as a label and form our MSE cost function calculate the difference between the predicted velocity and this data.
Acceleration (2 x 1000) matrix representing the acceleration. Not applicable our research, but could potentially be utilized for future research.
We will not experiment on the entire dataset of 30 motions due computational unfeasibility. Instead, select representative set of motions with (8 5 2 1000) samples in total. The experiments are mainly conducted with this designated set, but since the set is chosen to be representative, we expect the results to generalize to other motions as well.
Velocity normalization. Moreover, for some experiments, we opt to normalize the velocity values, such as in Figure 8, to avoid large cost values. This can cause loss of generality, since the policy actions are now restricted to the direction of the action vector, and will not try replicate its size. The size of this arrow might be important in scenarios where parts of the motion need be carried out at different pace. PLYDS can handle the dataset with or without velocity normalization. The dataset is also referenced and provided as part of our reproducibility efforts.
Real-world collected dataset. We collect data by teleoperating the Kinova Gen3 Lite Arm. Tele- operation involves employing human agents to operate a robotic device or system, recording their actions as expert demonstrations. The teleoperated actions are then recorded and utilized as training data for algorithmic learning. We have two options for teleoperation. First, a human expert can manually control the robot using joystick keyboard. This process results in natural but non-smooth trajectories. The second approach employs the robot's internal control systems and
14

---

Figure 8: Speed profile for the sine motion, normalized vs. natural. When we normalize the speed, policies fail capture the difference in the speed vector's magnitude along the trajectories. PLYDS works with both normalized and regular velocities, but we mostly opt for normalized velocities for baselining and comparisons, especially when plotting the policy streamlines for visual verification.
Polonged Sine Demoratratios
Pic and Poat Demoedratees
-2
12
Figure 9: Plots of the dataset collected using Kinova Gen3 Lite and teleoperation. The robot is operated complete the following trajectories multiple times, while the position and velocity data are recorded in real-time. Expert's 's demonstrations can also come from robot's 's low-level controllers, which leads faster data gathering process. We tried keep the scale of these trajectories aligned with the handwriting dataset to achieve consistency.
trajectory planning systems to perform as expert and execute some patterns. This leads to smoother data collection process with higher reliability. Please keep in mind that planners only connect a series of few points, with no guarantee of stability, and are time-dependent. Consequently, the role of policies generated by PLYDS will not obsolete.
We gather an open-source dataset holding three distinct motions: the prolonged sine, root parabola, and pick-and-plan (Figure 9). Each motion is represented by 50 demonstrations in 3-dimensional world. Each demonstration contains state (position) vector (3 1000) and corresponding action (velocity) vector (3 x 1000). Note that orientation is not recorded in the dataset, as we assume the robot's gripper will always face downwards. There will not be any loss of generality because controlling the orientation with PLYDS can be done in parallel, and in the exact same way as the end-effector's position. The dataset is provided as a part of our reproducibility efforts in Section 7.
B.2 SDP Optimization
We primarily use the commercially available MOSEK [42] optimization software that provides solutions for numerous types of optimization issues, including nonlinear semidefinite programming. The ficxibility and birf-peerormance capabilities of MOSEK make it ideal for challenging opti- mization tasks in both commercial and academic settings. We currently use the MOSEK under an academic license, which can be obtained free of charge with an academic domain email. SCS [48] is another solver specifically designed for solving semidefinite complementarity problems, which include nonlinear SDP special case. SCS employs augmented Lagrangian method combined
15

---

with the Burmeiter function to handle the complementarity conditions in the SDP. At this time, we do not have any solid comparison between the efficiency of these solvers for our setup, but commercial software products often perform more efficient than open-source products.
We also use SeiPy [49], an open-source scientific computing library for Python that has many modules for numerical optimization. SciPy handle wide range of optimization problems, including nonlinear programming with semidefinite constraints, even though it may not provide specialized solvers for nonlinear SDP. Our software still supports SciPy; however, it is not as efficient in solving nonlinear SDP problems as MOSEK and SCS.
B.3 Hyperparameters and architecture.
We provide a summary of parameters related to each baseline we used in the paper. Note that we accelerate the computation of GAIL, BC, and SDS-EF with NVIDIA GeForce RTX 3060 GPU, but SEDS, LPV-DS, and PLYDS use only Core-i7 Gen8 CPU for optimization.
PLYDS. For our experiments, we primarily utilize parameters a = 3 and = 1, which have proven effective in most cases. However, we also explore higher degrees to cover a broader range of settings. Additionally, to strike balance between stability and accuracy, we occasionally adjust the tolerance level from 10-1 to 10-9. This allows necessary. For the Lyapunov candidate maintain non-zero gradient beyond the quadratic form, we opt for only square elements in the LPF basis vector. This ensures stability and reliability in our system. Although it is possible to enforce a positive Hessian for the Lyapunov function, it incurs additional computational time while further limiting flexibility in stability conditions. More information about the parameters and architecture of PLYDS can be found on our GitHub repository: haabbb/ow/mmiiayaaaaaassaaiiiiii
trade off precision for stability when
The discriminator network takes as input the state-action pairs observations generated by the policy network and expert demonstrations. Hidden Layers: The network may consist of two or three hidden layers, each with 256 or 512 units. Activation Function: Rectified Linear Unit (ReLU) activation function is commonly used between the layers. Output: The discriminator produces a single output value, representing the probability of the input being from the expert or the generated policy. Hyperparameters:: Learning Rate: 0.0001, Number of Discriminator Updates per Generator Update: 1 or 2, Discount Factor (for reinforcement learning algorithm): 0.99, Batch Size: 64 or 128, Number of Training Iterations: 1000. We use the imitation package for GAIL's implementation: www.aaaaweeooo.oooee/////aaiiiie
GAIL..
The behavioral cloning network takes the state action pairs as input. Hidden Layers: The net- work may have two three hidden layers, each consisting of 128 256 units. Activation Function: Rectified Linear Unit (ReLU), Output: The output layer of the network corresponds the action space dimensionality, producing the predicted action. //yperparameters: Learning Rate: 0.0001, Number of Training Iterations: 5000, Batch Size: 64, Regularization Strength (L2 regularization): 0.001, Opti- mizer: Adam, Loss Function: Mean Squared Error (MSE). Same as with GAIL,, we use the imitation html.
BC.
package to access BC's implementation: aniatto.rrattedoosssc//////ttttr
SEDS. Takes cosition-velociny pairs input (or state-action pairs in general). Number of Gaussian Components: Typically ranges from 3 to 10, depending on the complexity of the motion being learned. Gaussian Mixture Model (GMM) Parameters: Covariance Type: Diagonal covariance is commonly used for efficiency and simplicity, Regularization Weight: Often set to a small value, such as 1e-6, to avoid singularities and overfitting, Maximum Number of Iterations: 100 iterations. Convergence Tolerance: le-6 1e-7. We have implemented SEDS in Python using the SeiPy optimization library, and the original MATLAB code not used in our comparisons to remain consistent with other baselines, particularly in terms of computational time. Our implementaaion of SEDS can be found on github. github.com/amina/ee http:///aaaaaaaaaaaaa-eeei-eeiia imilation-policy.
16

---

LPV-DS. We mainly use the original implementation of LPV-DS available on: githuh..oo/nbigggerrddd---e- and developed in Matlab. The parameters of the method re- main the same as the original repository, and we only change the number of demonstrations required for comparison purposes.
SDS-EF. We also provide an implementation of this baseline on our GitHub repository. The coupling layers set with the following parameters: base network 'rffn', activation func- tion = "elu', and sigma=0.45. The main architecture uses 10 blocks and hidden layers' size are set to 200. All these parameters are the same as SDS-EF's original implementation on gettubbcoommmaaaaatuudddddaaiiii
flows, but we omit the preprocessing step found in the origi- nal implementation, to be able to fairly, and effectively compare the results to other baselines.
C Supplemental Results
We present a comprehensive set additional experiments aimed at putting the proposed framework test from a variety of angles including access to fewer demonstrations, demonstrating the variety of LPFs, more baseline policy rollouts, additive noise, and lastly, we conduct ablation study by removing the stability guarantee, and present computation times in comparison with the baselines.
C.1 Baseline policy rollouts
In Figure 10, we plot policy rollouts optimized with PLYDS in comparison to the baselines: SEDS, GAIL,, and SDS-EF. We extend these results to LPV-DS Figure 11. All the portrayed policies are optimized on a set motions from the handwriting dataset, namely, G-Shaped, Angle, C-Shaped, and P-Shaped demonstrations. The key takeaway imitation learning methods for unknown areas of state space, and inaccuracies visible across the baselines. Hence, the same patterns the plots in the main text continue emerge.
pattern of instability among neural based
17

---

Pallcy Dipert start Eirpet
GOG
10
PLYDS
SEDS
GAAIL
-15
X1
Figure 10: Policy rollout for Angle, C-Shaped, G-Shaped, and P-Shaped demonstrations in hand- writing dataset. PLYDS is visually compared to the baselines in terms of reproduction accuracy and global stability. Note the inaccuracies and unstable reproductions in other baselines.
Cor
S
XI
Figure 11: Additional rollouts generated with the LPV-DS's source code in Matlab. These plots serve to complement the results acquired Figure 10.
C.2 Ablation study: stability vs. accuracy
When working with DS policies, there is a dilemma known as the stability-accuracy trade-off [22]. This means that a balance must be struck between the reliability and robustness of generated policies to guarantee global convergence to the target (referred to stability). and minimizing errors obtain precise solutions (referred to as accuracy). It is important to find compromise between these two
18

---

11
MAy
25
Unaaabe
23
Figure 12: Policy rollouts under utilizing PLYDS, both with and without the imposition of stability constraints, reveal significant difference policy behavior. The absence of enforced stability constraints, combined with the utilization of a complex polynomial and a preference for accuracy over stability embedded in the tolerance parameter, results system instability.
fapredksed
0.75
Esser Cempnetratirn
fium
-077
Deptee
Degree:
Degrve:
Figure 13: We delve into PLYDS policy rollouts, each employing distinct polynomial complexities. As evident from the plot, the use of increasingly complex polynomials results in more intricate trajectories that better mimic the expert's behavior, but generate more complex trajectories further from the demonstration data. This heightened complexity poses challenge in ensuring and validating the stability of the system.
factors, more stable algorithms may not be as accurate, while more accurate algorithms may be sensitive to instabilities.
The higher the degree of the polynomial, the more accurate imitation of expert behavior. Hence, in theory, any nonlinearity may be approximated by our DS formulation. However, practice, we introduce regularization and tolerance parameters in the code. The tolerance can be used choose in the favor of accuracy or stability (see Figure 12). Another way to balance this equation start with lower-degree polynomials, and increase the policy's complexity when the accuracy insufficient. Figure 13 serves as an illustration for this process.
C.3 Complexity of Lyapunov functions
Figure demonstrates the complexity of the Lyapunov function affects trajectory planning in the state space in various ways, such as optimization efficiency, trajectory smoothness, obstacle avoidance, robustness perturbations, and planning accuracy. If the Lyapunov function is more complex, may increase computational costs but in turn result in more complex and nonlinear trajectories. On the other hand, simpler Lyapunov functions may offer faster computations, smoother trajectories, and satisfactory planning accuracy, but they may not be adaptable to complex expert demonstrations.
When deciding on the complexity the Lyapunov function, necessary to consider the feasibility of computations, and the desired smoothness and accuracy of planned trajectories and always start with the most simple: quadratic distance function. Figure 14 illustrates this by showing both stable and unstable Lyapunov possibilities gauged across various LPF complexities.
19

---

Unasable
Subble
Degree
Degroe
Degroe
Figure 14: During our policy optimization, we obtain LPF samples such as the ones depicted above. Even though variation complexity notably affects the computation time, employing more complex Lyapunov functions (polynomials of higher degrees) appears necessary to achieve stable and precise policies some cases. Currently, we manually determine the complexity of the Lyapunov function and shift to higher complexities only if the optimization fails to deliver satisfactory results.
Refferecce Trajectoriess
inche 4.0)
tnctic 2.8)
incaue L8]
tncdec 0.5) (dda.l)
LegitSS)
X
Sine
Worm
D0line
Motion Shapes
Figure 15: Performance of PLYDS in the face of uniform additive noise (a) and sample of a noisy trajectory with noise-level set to 2 (b). Noise levels are in centimeters, therefore, a noise-level means each reading could be deviated from its true value by #4cm.
C.4 Performance with additive noise
Noise in imitation learning significantly impacts the learning process and resulting policies. Excessive noise levels can destabilize the algorithm, preventing it from converging to optimal policy. To assess PLYDS performance while exposed to noisy measurements, we apply uniform additive noise, distributing samples evenly across specified interval. We vary the size of this interval, expanding symmetrically around zero positions, while also accounting for its effect on velocities within the expert dataset. The results in Figure 15 demonstrate moderate level of noise-robustness that can be further improved in future studies. Noisy data also increases the error bands, leading to increased uncertainty in the outcome of policy optimization.
C.5 Computation times
We performed all experiments on machine equipped with Core i7 8th Gen CPU, an NVIDIA GeForce RTX 3060 GPU, and 32 GB DDR2 RAM. Among the methods included in experiments, GAIL, BC, and SDS-EF utilize the GPU expedite neural network computations. On the other hand,
20

---

400
PLYDS [dkg B)
15
PLYCS
PLYCS
PLYDS 100g G)
PLYDS
PLYDS (00g
Comput Contssaaioo Time (veconda))
-0-
PL/DS
300
COAL
230
200
alion
130
100
Sinc Worm Angli Motion Shapes
OBLinc
Sine worm Angle
DBUne
Motion Shapes
(b)
Figure 16: Total computation times averaged over 20 trials for PLYDS compared to other baselines (a) and with different dataset sizes (b). noteworthy that GAIL and BC are utilizing a GPU to accelerate their processing power.
FOTE krop (im
04
-.0
Mal Shipeo
13
(a)
(b)
(c)
Figure 17: Accuracy (MSE) and policy rollouts for vectorized (left) vs. scalar (right) Lyapunov function. For policy rollout, we picked the N N-Shaped motion, where the non-vectorized Lyapunov functions results in visible reduction in reproduction accuracy. However, the accuracy comparison shows that the extent of improvements caused by vector Lyapunov functions (marked by vec.) depend on the shape each motion, and may not be verified visually all motions.
PLYDS, SEDS, and LPV-DS solely rely the CPU for policy optimization. Despite this variance in computational resources, we provide comparison of computation times in Figure 16.
C.6
Ablation study: vector Lyapunov functions
Vector Lyapunov functions [37], extend the concept of Lyapunov stability to systems where state variables are represented as vectors rather than scalars. This extension particularly valuable when dealing with interconnected or multidimensional systems. Instead of relying on scalar function, our approach utilizes vector-valued Lyapunov function assign vector to each point in the state space. The properties of this vector-valued function are leveraged to analyze the stability and convergence behavior of system trajectories.
We employ this technique known to enhance the flexibility of the optimization process, as discussed [50]. Moreover, our observations indicate accuracy improvements over non-vectorized version for certain motion scenarios. Figure 17, we present the results of an ablation study focusing on vector Lyapunov functions, highlighting their essential role small yet critical component of our method. It is also known that Vectorizing the Lyapunov function yields higher flexibility during optimization, can potentially lower the computational cost for higher order polynomials.
21


Reviews:
[{"review_id": "co2WVYmlJ2", "reviewer": "Reviewer_ZV92", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "This paper proposes a method to learn stabilizing policy for polynomial dynamical systems using expert demonstrations. It parametrizes the policy via the P matrix and parametrizes the Lyapunov certificate via the Q matrix, and they jointly learn P and Q matrices via joint semi-definite programming (SDP). The optimization problem is solved via a classical SDP solver. They provide theoretical guarantees for stability (provided that the optimization is feasible). They conduct one simulation on handwriting and one hardware experiment for manipulation tasks. They outperform baselines (SEDS, BC, and GAIL) and get better performance for out-of-distribution samples.", "summary_of_recommendation": "Overall, I think the experiment has good initial results, but the paper's idea is a bit incremental and it lacks convincing baselines such as neural CLFs and GAIL with prior information for the goal point. Overall my decision is weak reject.", "technical_quality": "Fair: Significant additional evidence required to support the claims.", "clarity_of_presentation": "Good: Minor revision needed to improve organization or clarity.", "potential_impact": "", "robotics_focus": "Sufficient demonstration on hardware", "quality_of_limitations": "Limitations are addressed clearly", "questions_for_rebuttal": "I hope the authors can address more for the \"Cons\" section mentioned above. Besides are some extra questions:\n\n1. Why the Polynomial Lyapunov Candidate has to be designed? I read through section 4.2 but was not very clear about it. Since the final Lyapunov function is just a sum of all the Lyapunov candidates, why cannot just learn a Q=\\sum\\limits_{i=1}^n Q_i? Will it not work empirically?\n2. Lines 142-143 mentioned \"enables the policy and LPF to be learned independently for each dimension of the state-space\" - why those LPF are designed for each dimension, and why are those dimensions in the state space independent? From Eq.(4) it is hard to tell those dimensions are independent, and the state $\\mathtt{x}$ shows up in $\\mathtt{b_{x,\\beta}}$ and every dimension in $v(\\mathtt{x}; \\mathtt{Q})$ has the component of $\\mathtt{b_{x,\\beta}}$, thus certainly each dimension in the state space is not independent.\n3. Are all the expert trajectories not overlapping each other? Because in the state space, it is often required that from any point (other than the equilibrium point), there should be only one trajectory coming in (out).", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "**Pros:**\n1. The paper is written in clear English, and most sections are easy to understand.\n2. The experiments are conducted on non-trivial simulation and hardware cases, and the results are impressive.\n\n**Cons:**\n1. Do we really need Lyapunov/stability in imitation: the Lyapunov function only ensures the trajectories can finally approach the target point. But it cannot help more in mimicking the expert demonstration (other than reaching the set point), and the \"so-called\" stability is not that important in this case. One thought to make the method proposed here useful is to consider the stability in the error space - the error for the current state and the corresponding reference state in the expert trajectory. But from the current setup, it is unclear how easy it is to adapt to this \"error space\" setting.\n2. Baselines: In the related work, the author mentioned that neural Lyapunov functions typically lack GAS guarantees. However, the neural Lyapunov methods empirically still work well in high-dimension systems [1,2]. To fully show the advantage of the method proposed in this paper, I suggest comparing the neural Lyapunov methods. In addition, the BC and GAIL don't utilize the goal information as the proposed method and SEDS do, thus, the comparison is a little unfair. I suggest comparing it with [3], which utilizes the additional prior information (the goal point information). A simple comparison would be to enforce an L2 loss for the last several states in the planned trajectories and integrate it with the loss function in Eq (4) in [3]. Also, in terms of neural network structure, a recurrent neural network might be preferable since it captures better the temporal information of the trajectory.\n\nReferences:\n1. Dawson, Charles, et al. \"Safe nonlinear control using robust neural lyapunov-barrier functions.\" Conference on Robot Learning. PMLR, 2022.\n2. Chang, Ya-Chien, and Sicun Gao. \"Stabilizing neural control using self-learned almost Lyapunov critics.\" 2021 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2021.\n3. Bhattacharyya, Raunak P., et al. \"Simulating emergent properties of human driving behavior using multi-agent reward augmented imitation learning.\" 2019 International Conference on Robotics and Automation (ICRA). IEEE, 2019.", "extra_scores": {"impact": "3: The work is incremental and unlikely to have much impact, even though it may be technically correct and well executed."}, "extra_text": {}}, {"review_id": "Utuq-xjV00", "reviewer": "Reviewer_x3iz", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "This paper presents a method for learning stable dynamical systems from demonstrations. To obtain hard guarantees on stability, the learned dynamics and Lyapunov stability certificate are represented with polynomials, and an optimization is solved to synthesize the dynamics and Lyapunov function jointly. The method is evaluated on the LASA dataset and on manipulation tasks, outperforming baselines in stability and demonstration tracking.", "summary_of_recommendation": "The paper is overall well-written but needs some clarification on the method, expansion of related work, and comparison with some more modern baselines.", "technical_quality": "Good: Generally convincing, but some additional experiments or analyses needed.", "clarity_of_presentation": "Good: Minor revision needed to improve organization or clarity.", "potential_impact": "", "robotics_focus": "Sufficient demonstration on hardware", "quality_of_limitations": "Additional details required", "questions_for_rebuttal": "Please see review for specific questions.", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Strengths: The paper is well-written, has extensive experiments, and as a contribution in stable learning from demonstration is of relevance to the robot learning community.\n\nWeaknesses and concerns:\n- There is a gap in related work and a lack of comparison to modern approaches for learning stable systems with guarantees.\n    - The idea of using leveraging polynomials to learn and certify stability of dynamics is not a new one. For instance, consider the work [1] on learning contracting dynamical systems along demonstration trajectories with local certificates of convergence, or [2] for learning Lyapunov stable polynomial latent dynamics in the context of control from pixels.\n    - Recent methods leveraging ideas from Riemannian motion policies [3, 4] which learn diffeomorphisms from a known stable policy to a latent space with dynamics that are trained to resemble the demonstrations.\n    - I would say the claim that neural Lyapunov methods “typically lack GAS since they softly enforce the Lyapunov stability conditions” in the related work is a bit misleading. Despite using NNs for the dynamics and Lyapunov function, [5] ([33] in the paper) does provide hard stability guarantees (for the learned dynamics) by using counterexample-guided synthesis, and verification using mixed integer programming. This provides hard guarantees at the cost of an expensive verification procedure. On the other end of the spectrum, there are methods leveraging the Lipschitz constants of the Lyapunov/contraction conditions of the learned certificate function and using these to certify violations of the conditions cannot occur over some compact set [6, 7]; these approaches are more conservative in verification but are more scalable than mixed integer verifiers.\n- I am confused as to how the optimization in (7) is a (convex) SDP. If one searches for the dynamics (parameterized by P) and the Lyapunov function (parameterized by Q), there should be bilinear terms in the coefficient matching constraint (5), which renders the problem non-convex. There are approaches for solving this iteratively through alternating between fixing P and optimizing Q and vice versa (like in [2]), but it seems this should be a non-convex problem at its core.\n- The results would be more convincing if the method is compared with the existing techniques in learning stable polynomial systems from demonstrations (e.g., [1] or [4]).\n- The global stability guarantees coming from the SOS program are at odds with the limitations of practical robotic systems, which have bounded domains and control limits. This can compromise the guarantees of the learned behavior. For instance, since the Kinova arm in the experiments has joint limits, there could be the case where the arm is perturbed to some point near the the joint limits where the learned dynamical system drives the arm to exceed the joint limits in an unsafe fashion.\n- It would be good to see an ablation on using vector Lyapunov vs. scalar Lyapunov functions. I am unclear on the benefit of having separate Lyapunov functions per dimension and then enforcing that they elementwise all satisfy the standard Lyapunov conditions. It seems that if you can find Q_1, …, Q_n which all satisfy the Lyapunov conditions independently, you should be able to find a single Q for a scalar Lyapunov function (maybe with a polynomial basis of slightly higher degree). I think the idea of using vector Lyapunov functions is interesting, as they admit weaker stability conditions (e.g., component-wise you do not need to be globally negative semidefinite, see [8]), but it seems their advantages are not yet being fully exploited in this work.\n\n[1] El Khadir, Varley, Sindhwani. Teleoperator Imitation with Continuous-time Safety, RSS’19.\n\n[2] Chou and Tedrake. Synthesizing Stable Reduced-Order Visuomotor Policies for Nonlinear Systems via Sums-of-Squares Optimization, CDC’23.\n\n[3] Asif Rana et al., Euclideanizing Flows: Diffeomorphic Reduction for Learning Stable Dynamical Systems, L4DC’20\n\n[4] Zhang, Beik-Mohammadi, and Rozo. Learning Riemannian Stable Dynamical Systems via Diffeomorphisms, CoRL’22.\n\n[5] Dai, Landry, Yang, Pavone, and Tedrake. Lyapunov-stable neural network control, RSS’21.\n\n[6] Chou, Ozay, and Berenson. Model Error Propagation via Learned Contraction Metrics for Safe Feedback Motion Planning of Unknown Systems, CDC’21.\n\n[7] Sun, Jha, and Fan. Learning Certified Control using Contraction Metric, CoRL’20.\n\n[8] Nersesov and Haddad. On the Stability and Control of Nonlinear Dynamical Systems via Vector Lyapunov Functions, TAC’06.", "extra_scores": {"impact": "3: The work is incremental and unlikely to have much impact, even though it may be technically correct and well executed."}, "extra_text": {}}, {"review_id": "Eifkx7WMl5", "reviewer": "Reviewer_fdCw", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "The authors propose the use of stable polynomial dynamical systems to encode expert demonstrations in the context of learning from demonstrations.\n\nTo guarantee the stability and accuracy of the polynomial dynamical systems, the authors propose an optimization method based in sum-of-squares optimization. \n\nThe paper includes a thorough comparison with state-of-the-art methods, and the experiments demonstrate that this method can outperform the state-of-the-art in specific scenarios while also being efficient in terms of sample usage.\n", "summary_of_recommendation": "While polynomial dynamical systems have been extensively studied in the control and systems literature, their application in the context of motion planning has not been explored until now.\n\nThis paper fills that gap by connecting polynomial dynamical systems with motion planning, thereby paving the way for further research. By incorporating the substantial body of literature on sum-of-squares (SOS) methods and polynomial dynamical systems into motion planning, this paper makes a valuable contribution to the field of learning from demonstration (LfD) using dynamical systems.", "technical_quality": "Good: Generally convincing, but some additional experiments or analyses needed.", "clarity_of_presentation": "Very Good:  Well organized and clear.", "potential_impact": "", "robotics_focus": "Sufficient demonstration on hardware", "quality_of_limitations": "Additional details required", "questions_for_rebuttal": "- ε in (2) is not defined \n\n- It is said that (7) is a semi-definite programming with a nonlinear cost function, but constraint (c) of (7) could induce bilinear/equality constraints. When this happens, (7) is not more an SDP in its constraints nor a convex problem.\n\n- No reference is cited for the SCP solver\n\n- It is claimed that higher degrees than the ones used in the authors' experiment are often unnecessary in normal tasks, but in Appendix B.1 is also said that the method was not experimented on the entire dataset due to computational unfeasibility. Please clarify.\n\n- In the Limitations subsections, it said that the Lasserre hierarchy and SOS extensions presented in [43] could search in a broader class of LPF candidates and tackle the issue of finding Lyapunov functions. However, the methods presented in [43] assume that f is given, while in the paper, a simultaneous search needs to be also carried out to find $f$. Thus, in the reviewer's opinion, this does not directly solve the issue. Nevertheless, the authors could check the methods presented in \"A Survey of Recent Scalability Improvements for Semidefinite Programming with Applications in Machine Learning, Control, and Robotics\" which aims to mitigate the scalability of SOS methods in general.\n\n- The authors write that \"A limitation of SOS is that the set of non-negative polynomials is larger than the ones expressed as SOS [42]. Though rare in motion planning tasks, this implies that finding a Lyapunov candidate could be difficult\". This implication is not entirely true, as Ahmadi and Parrilo proved in \"Converse results on existence of sum of squares Lyapunov functions\" that the existence of a polynomial Lyapunov function of a given degree implies the existence of a polynomial Lyapunov function $V$ of (possibly higher degree) that satisfies\nthe sos conditions $V$ is SOS and $-dV/dt$ is SOS. Nevertheless, the authors are simultaneously searching for \\hat{f} and $V$, so the limitation presented in [42] is relevant here. But I think the authors should clarify that the result of [42] implies that simultaneously finding the polynomial system and the Lyapunov function could be difficult. \n\n- It is known that the scalability of vanilla SOS methods is very challenging since the size of the matrices involved in the problem can quickly grow with the degree. This should be more clearly stated in the Limitation subsection.  \n", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper has several strengths, including comprehensive comparisons with various methods and demonstrations on a robotic platform. Additionally, the authors provide the implementation code, which facilitates result reproducibility and enhances the paper's impact.\n\nHowever, there are a few drawbacks to the method presented in the paper. For instance, the scalability of the method is limited, as it was not tested on the entire LASA dataset due to computational unfeasibility. Furthermore, certain aspects of the paper could benefit from further clarification.\n", "extra_scores": {"impact": "4: The work contains interesting new ideas that will potentially have major impact in robotics or machine learning."}, "extra_text": {}}, {"review_id": "-IZS31hdL2z", "reviewer": "Reviewer_gRUM", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "This paper offers a solution to the stability vs. accuracy dilemma that is the crux of learning reactive motion policies as stable dynamical systems from demonstrations. In this work, this is achieved by formulating both the DS motion policy and the corresponding Lyapunov function (used to enforce stability) as multidimensional polynomial functions. An issue of Lyapunov-based DS learning approaches is the need to estimate both the policy and the Lyapunov functions. In this paper, the parameters of both of these polynomial functions are estimated jointly using SOS by generating a single polynomial representation constructed from a polynomial coefficient matching strategy. This presented approach, named PLYDS, is shown to outperform the first baseline stable DS learning method from 2012 (SEDS) and unstable learning approaches like BC and GAIL, in terms of reproduction accuracy while ensuring global asymptotic stability.", "summary_of_recommendation": "The paper presents a novel polynomial dynamical system formulation and efficient learning scheme to tackle the stability vs. accuracy dilemma. The approach is sound, however, the comparisons are insufficient to support the claims. If the authors can address the weaknesses I mention above and particularly show improved performance over any (at least one) of the state-of-the-art approaches that outperform SEDS I would be willing to increase my score. \n\n**Post rebuttal**\n\nThe authors have addressed my major concerns regarding comparisons to state-of-the-art. I recommend they highlight these comparisons in the main text. Also, even though it was not addressed in the rebuttal it would be great if the authors could discuss the scalability of their approach to higher dimensional systems and manifolds such as SE(3) and joint space.", "technical_quality": "Good: Generally convincing, but some additional experiments or analyses needed.", "clarity_of_presentation": "Very Good:  Well organized and clear.", "potential_impact": "", "robotics_focus": "Sufficient demonstration on hardware", "quality_of_limitations": "Additional details required", "questions_for_rebuttal": "Please address weaknesses listed above.", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "**Strengths**\n- The paper is well-written and enjoyable to read. \n- I went through all the proofs in the appendix and the math seems sound and refreshingly very well explained. \n- This approach can be seen as is a very interesting application of \"polynomial regression\", which is having a resurgence thanks to recent comparisons to NN. \n- The transformation of P and Q to G is a very clever trick to solve the problem of learning the parameters of each of the polynomial functions as it allows them to use standard SOS to estimate them.\n- I did not have time to test the code, but it looks very well written too.\n- Solid solution to a well-known problem.\n\n**Weaknesses**\n1. **Motivation/Justification:** Although I appreciate this polynomial formulation and the fact that it allows efficient parameter estimation, I am yet to understand why another SEDS-improvement approach is necessary that only solves for the stability vs. accuracy dilemma. As mentioned by the authors, there are many improvements to SEDS, like the \n   - i) the CLF learning approach proposed in 2015  https://www.sciencedirect.com/science/article/abs/pii/S0921889014000372\n   - ii) the LPV-DS approach that decouples the GMM parameter estimation to the linear DS estimation stages and uses a P-shaped Lyapunov function already shows superior performance for the type of motions that are being demonstrated in this work (highly non-linear, single-attractor systems) https://proceedings.mlr.press/v87/figueroa18a.html\n   - iii) the Euclideanizing-Flow approach that learns a diffeomorphic transformation using normalizing flows https://proceedings.mlr.press/v120/rana20a.html\n\n   I only mention these 3 as they are the ones I know that have implementations/code available offline for the authors to compare.  If the PLYDS approach gives the same performances as these improved (more recent) solutions, then what is the contribution from this new formulation? In general, the majority of the state-of-the-art approaches can already achieve similar reproduction accuracy on the LASA dataset and can handle highly nonlinear trajectories. They all, however, suffer from either computational complexity (during learning and execution) as well as hyper-parameter tuning. So, can the authors clarify or spell out what is the improvement of PLYDS not against SEDS but against the more recent state-of-the-art that already outperform SEDS? This leads me to my next point.\n\n2. **Insufficient Evaluation:**  Using SEDS as a baseline + knowingly unstable policies (BC and GAIL) is not a valid comparison in my opinion. The work presented in this paper is offering an improvement on SEDS, yet many improvements of SEDS exist. It's fine to have SEDS as the lowest baseline approach which we know will not perform well, but to corroborate your claims you should at least compare to one of the improved DS learning approaches to have a more fair performance evaluation. Browsing through your code I saw that SEDS is implemented from scratch in Python. The authors could consider implementing this approach http://proceedings.mlr.press/v78/ravichandar17a/ravichandar17a.pdf which is exactly the same SEDS formulation but with partial contraction constraints instead of Lyapunov. Furthermore, there is no discussion about the computational time in the main text, I see that this is discussed in the appendix. If possible, I suggest to move some of this relevant information to the main text as it's one of important improvement. Can you clarify if for the SEDS computation time is considering the model selection part or is it only timing the SEDS learning algorithm on a predefined number of Gaussians?  This leads to the next point.\n\n3. **Hyper-parameter Tuning:** I feel like the fact that the complexity of the polynomials for both the motion policy and the Lyapunov function are manually defined is not fully addressed in the main text. As you already show in the appendix. the degree of the polynomial affects both reproduction accuracy and stability so a proper set of degrees for each polynomial function must be found, otherwise the claim that PLYDS will outperform SEDS does not hold universally. However, this is not such a pressing issue as very few of the DS learning approaches (with the exception of https://proceedings.mlr.press/v87/figueroa18a.html) requires model selection or hyper-parameter grid search to fund such optimal parameters. Nevertheless, it would be great if the authors discussed or provided some rule-of-thumb for selecting these hyper-parameters or an approach to do such selection automatically or in a data-driven way.\n\n3. **Scalability to higher dimensions:** The experiments in this paper are limited to 2D-3D position trajectories of the robot's end-effector. As I mentioned before, learning stable DS in 2D-3D from highly nonlinear trajectories has already been solved in a myriad of approaches. Yet, apart from the computational limitations many formulations also lack the ability to scale to higher dimensional complex spaces, like $SE(3)$ full end-effector pose or in joint space $q\\in R^N$, like RMPs https://arxiv.org/abs/1801.02854. Hence, a more relevant improvement would be to show that this polynomial formulation allows to learn DS in such high-dimensional spaces. Is it possible? How would this affect the required degree of polynomials and computation efficiency, if at all?", "extra_scores": {"impact": "3: The work is incremental and unlikely to have much impact, even though it may be technically correct and well executed."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper tackles the problem of imitation learning under stability guarantees by learning a parametric polynomial planning function. The approach learns the polynomial coefficients jointly with a Lyapunov candidate to ensure stability.

All reviewers appreciated the clarity and the extensive experimental evaluation of the paper, particularly after the author's response that clarified the remaining doubts and presented additional comparisons.

the authors' effort in the rebuttal was appreciated and led the reviewers to increase their scores.

For the camera-ready version, I suggest the authors improve the related work's sections, as many important citations are missing, particularly on the diffeomorphisms-based approach to stability and stability on the Riemannian manifolds (only citing the Euclideanizing Flows is not enough, as many other works in the area exist). Furthermore, it would be nice to include an in-depth discussion on the scalability of the method, particularly w.r.t. the diffeomorphism-based approaches.
Decision: Accept (Poster)