On the Sample Complexity of Stability Constrained Imitation Learning

Introduction

Imitation Learning (IL) techniques use demonstrations of desired behavior, provided by an expert, to train a policy. IL offers many appealing advantages: it is often more sample-efficient than reinforcement learning , and can lead to policies that are more computationally efficient to evaluate online than optimization-based experts. Indeed, there is a rich body of work demonstrating the advantages of IL-based methods in a range of applications including video-game playing , humanoid robotics , and self-driving cars . Safe IL further seeks to provide guarantees on the stability or safety properties of policies produced by IL. Methods drawing on tools from Bayesian deep learning , PAC-Bayes , stability regularization , or robust control , are able to provide varying levels of guarantees in the context of IL.

However, when applied to continuous control problems, little to no insight is given into how the underlying stability properties of the expert policy affect the sample-complexity of the resulting IL task. In this paper, we address this gap and answer the question: what makes an expert policy easy to learn? Our main insight is that when an expert policy satisfies a suitable quantitative notion of robust incremental stability, i.e., when pairs of system trajectories under the expert policy robustly converge towards each other, and when learned policies are also constrained to satisfy this property, then IL can be made provably efficient. We formalize this insight through the notion of incremental gain stability constrained IL algorithms, and in doing so, quantify and generalize previous observations of efficient and robust learning subject to contraction based stability constraints.

There exist a rich body of work examining the interplay between stability theory and learning dynamical systems/policies satisfying stability/safety properties from demonstrations.

Nonlinear stability and learning from demonstrations: Our work applies tools from nonlinear stability theory to analyze the sample complexity of IL algorithms. Concepts from nonlinear stability theory, such as Lyapunov stability or contraction theory , have also been successfully applied to learn autonomous nonlinear dynamical systems satisfying desirable properties such as (incremental) stability or controllability. As demonstrated empirically in , using such stability-based regularizers to trim the hypothesis space results in more data-efficient and robust learning algorithms. However, no quantitative sample-complexity bounds are provided.

To provide fine-grained insights into the relationship between system stability and sample-complexity, we first define and analyze the notion of incremental gain stability (IGS) for a nonlinear dynamical system. IGS provides a quantitative measure of robust convergence between system trajectories, that in our context strictly expands on the guarantees provided by contraction theory by allowing for a graceful degradation away from exponential convergence rates.

We then propose and analyze the sample-complexity properties of IGS-constrained imitation learning algorithms, and show that the graceful degradation in stability translates into a corresponding degradation of generalization bounds by linking nonlinear stability and statistical learning theory. In particular, we show that when imitating an IGS expert policy, IGS-constrained behavior cloning requires $m\gtrsim q\cdot T^{2a(1-1/a^{2})}\cdot\varepsilon^{-2a}$ trajectories to achieve imitation loss bounded by $\varepsilon$, where $T$ is the task horizon, $q$ is the effective number of parameters of the function class for the learned policy, and $a\in[1,\infty)$ is an IGS parameter determined by the expert policy. We show that $a=1$ for contracting systems, leading to task-horizon independent bounds scaling as $m\gtrsim{q}/{\varepsilon^{2}}$. Furthermore, we construct a simple family of systems where the IGS parameter $a$ satisfies $a=1+p$ for $p\in(0,\infty)$. This yields sample-complexity that scales as $m\gtrsim q\cdot T^{\nicefrac{{2p(2+p)}}{{1+p}}}\cdot\varepsilon^{-2(1+p)}$, which makes clear that an increase/decrease in $p$ yields a corresponding increase/decrease in sample-complexity.

Motivated by the empirical success and widespread adoption of DAgger and DAgger-like algorithms, we also extend our analysis to an IGS-constrained DAgger-like algorithm. We show that this algorithm enjoys comparable stability dependent sample-complexity guarantees, requiring $m\gtrsim q\cdot T^{2a^{2}\left(1-1/a\right)\left(1+1/a+3/2a^{2}\right)}\cdot\varepsilon^{-2a^{2}}$ trajectories to achieve $\varepsilon$-bounded imitation loss, again recovering time-independent bounds for contracting systems that gracefully degrade when applied to our family of systems satisfying $a=1+p$.

Together, our results are the first to delineating a class of systems where the sample-complexity bounds for imitation learning scale sublinearly in the task-horizon $T$, and do so without requiring (strong) convexity of the loss function in the policy parameters. We conclude by demonstrating the validity of our theoretical results on (a) our simple family of nonlinear systems for which the underlying IGS properties can be quantitatively tuned, and (b) a high-dimensional nonlinear quadrupedal robotic system. Empirically, we find that the sample-complexity scaling predicted by the underlying stability properties of the expert policy are indeed observed in practice.

Problem Statement

We consider the following discrete-time dynamical system:

Incremental Gain Stability

The crux of our analysis relies on a property which we call incremental gain stability (IGS). Before formally defining IGS, we motivate the need for a quantitative characterization of convergence rates between system trajectories. A key quantity that repeatedly appears in our analysis is the following sum of trajectory discrepancy induced by policies $\pi_{1}$ and $\pi_{2}$:

We already saw this quantity appear naturally in (2.3). Furthermore, we will reduce analyzing the performance of behavior cloning and our DAgger-like algorithm to bounding the discrepancy (3.1) between trajectories induced by the expert policy $\pi_{\star}$ and a learned policy $\hat{\pi}$.

The simplest way to bound (3.1) is to use a discrete-time version of Grönwall’s inequality: if the map $f$ defining (2.1) is $L_{f}$-Lipschitz, in addition to the policies $\pi_{1}$ and $\pi_{2}$ being $L_{\pi}$-Lipschitz, then (assuming $L_{f}L_{\pi}\gg 1$) we can upper bound the discrepancy (3.1) by:

allows us to treat $\{\Delta_{\pi_{1},\pi_{2}}(x_{t}^{\pi_{1}}(\xi))\}_{t\geqslant 0}$ as an input signal, yielding

Here, $\Delta_{t}:=x_{t}(\xi_{1},\{u_{t}\}_{t\geqslant 0})-x_{t}(\xi_{2},\{0\}_{t\geqslant 0})$.

IGS quantitatively bounds the amplification of an input signal $\{u_{t}\}_{t\geqslant 0}$ (and differences in initial conditions $\xi_{1},\xi_{2}$) on the corresponding system trajectory discrepancies $\{\Delta_{t}\}_{t\geqslant 0}$. Note that a system that is incrementally gain stable is automatically $\delta$ISS. IGS also captures the phase transition that occurs in non-contracting systems about the unit circle. For example, when $\|\Delta_{t}\|\leqslant 1$ and $\|u_{t}\|\leqslant 1$ for all $t\geqslant 0$, inequality (3.3) reduces to $\sum_{t=0}^{T}\lVert\Delta_{t}\rVert_{X}^{a}\leqslant\zeta\lVert\xi_{1}-\xi_{2}\rVert_{X}^{\alpha_{0}}+\gamma\sum_{t=0}^{T-1}\lVert u_{t}\rVert_{U}.$ Finally, as IGS measures signal-to-signal ($\{u_{t}\}_{t\geqslant 0}\to\{\Delta_{t}\}_{t\geqslant 0}$) amplification, it is well suited to analyzing learning algorithms operating on system trajectories.

Then, $f$ is $(a,b,\Psi)$-incrementally-gain-stable with $\Psi=\left(\alpha_{0},\frac{\overline{\alpha}}{\underline{\alpha}\wedge\mathfrak{a}},\frac{\mathfrak{b}}{\underline{\alpha}\wedge\mathfrak{a}}\right)$.

Our first example of incremental gain stability is a contracting system .

Consider the dynamics $x_{t+1}=f(x_{t},u_{t})$. Suppose that $f$ is autonomously contracting, i.e., there exists a positive definite metric $M(x)$ and a scalar $\rho\in(0,1)$ such that:

Three concrete examples of autonomously contracting systems include:

$f(x,u)=\log(1+x^{2})+u$ with the metric $M(x)=2[1+\exp(-|x|)]^{-1}$, and

$f(x,u)=x-\eta\left[\nabla V(x)+u\right]$ where $V(x)$ is a twice differentiable potential function satisfying $\mu I\preccurlyeq\nabla^{2}V(x)\preccurlyeq LI$, and $0<\eta\leqslant 1/L$ .

Our next example illustrates a family of systems that degrade away from exponential rates.

Consider the scalar dynamics $x_{t+1}=x_{t}-\eta x_{t}\frac{|x_{t}|^{p}}{1+|x_{t}|^{p}}+\eta u_{t}$ for $p\in(0,\infty)$. Then as long as $0<\eta<\frac{4}{5+p}$, we have that $f$ is $(1+p,1,\Psi)$-IGS, with $\Psi=\left(1,\frac{2^{2+p}}{\eta},2^{2+p}\right)$

The system described in Proposition 3.5 behaves like a stable linear system when $|x_{t}|\geqslant 1$, and like a polynomial system when $|x_{t}|<1$ (hence $a=1+p$). This example highlights the need to be able to capture a phase-transition within our definitions and Lyapunov characterizations.

Algorithms and Theoretical Results

In this section we define and analyze IGS-constrained imitation learning algorithms. We begin by introducing our main assumption of dynamics and policy class regularity.

We assume that the dynamics $f$, policy class $\Pi$, expert $\pi_{\star}$, and initial condition distribution $\mathcal{D}$ satisfy:

The dynamics map $f$ satisfies $f(0,0)=0$.

The policy class $\Pi$ is convex and $\pi(0)=0$ for all $\pi\in\Pi$.

The distribution $\mathcal{D}$ over initial conditions satisfies $\lVert\xi\rVert_{2}\leqslant B_{0}$ a.s. for $\xi\sim\mathcal{D}$.

$\Delta_{\pi_{1},\pi_{2}}$ is $L_{\Delta}$-Lipschitz for all $\pi_{1},\pi_{2}\in\Pi$.

The constants $B_{0},L_{\Delta}\in[1,\infty)$.

Before turning to our main stability assumption we briefly remark on Assumption 4.1(b), which requires that the policy class $\Pi$ is convex. This assumption is stronger than is actually necessary. Instead, we could consider optimizing at epoch $k$ over a function class $\Pi_{k}$ defined recursively as $\Pi_{k}=\{\alpha\pi+(1-\alpha)\pi_{k-1}\mid\alpha\in,\,\pi\in\Pi,\,\pi_{k-1}\in\Pi_{k-1}\}$, with the base case $\Pi_{0}=\Pi$. Instead, we choose to make the assumption that $\Pi$ is convex to simplify the presentation, noting that using the recursive representation would yield the same sample complexity bounds.Our results are derived by bounding the Rademacher complexity of a particular function class, which in our setting is preserved under convex hulls. See Proposition E.4 in the appendix for more details.

We remark that we assume that the expert policy lies in our policy class, i.e., $\pi_{\star}\in\Pi$, to guarantee that zero imitation loss can be achieved in the limit of infinite data; it is straightforward to relax this assumption to $\Pi\cap\mathcal{S}(a,1,\Psi)\neq\emptyset$ and prove results with respect to the best stabilizing policy in class.

With these definitions and assumptions in place, we introduce IGS-Constrained Mixing Iterative Learning (CMILe) in Algorithm 1 and state our main theoretical results. CMILe draws upon and integrates ideas from Stochastic Mixing Iterative Learning (SMILe) , constrained policy optimization , and the IGS tools developed in Section 3. As in SMILe and DAgger, CMILe proceeds in epochs, beginning with data generated by the expert policy, and iteratively shifts towards a learned policy via updates of the form $\pi_{k+1}=(1-\alpha)\pi_{k}+\alpha\hat{\pi}_{k}$, where $\pi_{k}$ is the current data-generating policy, $\hat{\pi}_{k}$ is the policy learned using the most recently generated data, and $\alpha\in(0,1]$ is a mixing parameter. However, CMILe contains two key departures from traditional IL algorithms: (i) it constrains the learned policy at each epoch to remain appropriately close to the previous epoch’s data-generating policy (constraint (4.1b)), and (ii) all data-generating policies $\{\pi_{k}\}$ are constrained to induce IGS closed-loop systems (constraint (4.1c)). The latter constraint allows us to leverage the IGS machinery of Section 3 to analyze Algorithm 1.

In presenting our results, we specialize the policy class $\Pi$ to have the parametric form:

with $B_{\theta}\geqslant 1,$ and $\pi$ a fixed twice continuously differentiable map. As an example, neural networks with $q$ weights and twice continuously differentiable activation functions are captured by the policy class (4.2). We note that our results do not actually require a parameteric representation: as long as a particular policy class Rademacher complexity (defined in Appendix E) can be bounded, then our results apply. In what follows, we define the following constants:

We first analyze a single epoch version of Algorithm 1, which reduces to Behavior Cloning (BC) subject to interpolating the expert policy on the training data (constraint (4.1b)) and inducing an IGS closed-loop system (constraint (4.1c)).

Suppose that Assumption 4.1 and Assumption 4.2 hold. Set $\alpha=E=1$ in Algorithm 1. Suppose that $m$ satisfies:

With probability at least $1-e^{-q}$ over the randomness of Algorithm 1, we have that:

Theorem 4.3 shows that the imitation loss for IGS-constrained BC decays as $T^{1-1/a^{2}}\cdot(\tfrac{q}{m})^{\frac{1}{2a}}$ We discuss implications on sample-complexity after analyzing the general setting.

Next we analyze Algorithm 1 as stated, and show that if the mixing parameter $\alpha$ and number of episodes $E$ are chosen appropriately with respect to the IGS parameters of the underlying expert system, sample-complexity guarantees similar to those in the IGS-constrained BC setting can be obtained. As described above, the key to ensuring that guarantees can be bootstrapped across epochs is the combination of a trust-region constraint (4.1b) and IGS-stability constraints (4.1c) on the intermediate data-generating policies.

Suppose that Assumption 4.1 and Assumption 4.2 hold, and that:

Suppose further that for $k\in\{1,\dots,E-2\}$, we have:

that $E$ divides $m$, $E\geqslant\frac{1}{\alpha}\log\left(\frac{1}{\alpha}\right)$, and $\alpha\leqslant\min\left\{\frac{1}{2},\frac{1}{L_{\Delta}\gamma T^{1-1/a}}\right\}$. Then with probability at least $1-e^{-q}$ over the randomness of Algorithm 1, Algorithm 1 is feasible for all epochs, and:

Theorem 4.4 states that if the mixing parameter $\alpha$ and number of episodes $E$ are set according to the underlying IGS-stability parameters of the expert system then the imitation loss of the final policy $\pi_{E}$ scales as $T^{\left(1-\frac{1}{a}\right)\left(1+\frac{1}{a}+\frac{3}{2a^{2}}\right)}\cdot\left(\frac{q}{m}\right)^{\tfrac{1}{2a^{2}}}$.

By comparing the sample complexity bound for IGS-BS (Theorem 4.3) to the bound for IGS-CMILe (Theorem 4.4), we see that for fixed IGS parameters and number of trajectories $m$, the imitation error for IGS-BC is order-wise dominated by the IGS-CMILe imitation error. That is, while our current analysis does show the benefits of expert robustness via explicit dependence on the IGS parameters, it does not show the relative benefits of IGS-CMILe over IGS-BC, despite our experimental evidence suggesting otherwise (cf. Section 5). We leave a theoretical analysis showing the benefit of IGS-CMILe over IGS-BC to future work.

From the above discussion, we can delineate classes of systems for which imitation loss sample-complexity bounds are sublinear in the task horizon $T$. Specifically, we bound the number of trajectories $m$ needed to achieve $\varepsilon$-bounded imitation loss, ignoring logarithmic factors and problem constants except the horizon length $T$ and the IGS parameter $a$:

IGS-BS (Theorem 4.3) requires $m\gtrsim\varepsilon^{-2a}\cdot T^{2a(1-1/a^{2})}$ trajectories; this is sublinear in $T$ when $a\in[1,(1+\sqrt{17})/4)\approx[1,1.281)$.

IGS-CMILe (Theorem 4.4) requires $m\gtrsim\varepsilon^{-2a^{2}}\cdot T^{2a^{2}\left(1-\frac{1}{a}\right)\left(1+\frac{1}{a}+\frac{3}{2a^{2}}\right)}$ trajectories; this is sublinear in $T$ when $a\in[1,(3/2)^{1/3})\approx[1,1.144)$.

Finally, when a system is contracting, we have $a=1$ by Proposition 3.4, and hence the number of required trajectories $m$ for both IGS-BC and IGS-CMILe simplifies to $m\gtrsim\varepsilon^{-2}$.

The requirement on the constraint slack $c_{k}$ in (4.4) allows the constrained ERM problem (Algorithm 2) non-zero slack in matching the behavior of the previous policy (cf. (4.1b)). This is compatible with practical implementations of first order trust region policy optimization, where constraints are enforced via soft losses instead of as hard constraints.

1 Necessity of Stability Constraints

The empirical risk minimization algorithm (Algorithm 2) we consider in this work requires an IGS constraint on the learned policy (cf. (4.1c)). Here, we show the necessity of imposing this stability constraint in order to derive high probability sub-exponential in $T$ bounds on the imitation error. Consider the linear time-invariant system:

which contains $\pi_{\star}$. Note that Assumptions 4.1 and 4.2 hold for the closed-loop expert dynamics and policy class. The behavior cloning ERM problem (without stability constraints on $K$) is:

It is easy to check that on $\mathcal{E}_{m}$ (a constant probability event), $\hat{K}=-2e_{1}e_{1}^{\mathsf{T}}$ is an optimal solution of this ERM problem (that achieves zero training loss). Let $\hat{\pi}(x)=\hat{K}x$. Since $x_{t}^{\hat{\pi}}(e_{2})=2^{t}e_{2}$,

Therefore, if one removes the stability constraint (4.1c), then sub-exponential in $T$ imitation error bounds are impossible without more problem assumptions or other algorithmic modifications.

Experiments

In our experiments, we implement neural network training by combining the haiku NN library with optax in jax .

In order to implement Algorithm 1, a constrained ERM subproblem (Algorithm 2) over the policy class must be solved. Two elements make this subproblem practically challenging: (i) the trust-region constraint (4.1b), and (ii) the IGS-stability constraint (4.1c).

We implement the trust-region constraint (4.1b) by initializing the weights $\hat{\theta}_{k}$ parameterizing the policy $\hat{\pi}_{k}$ at those of the previous epoch’s parameters $\hat{\theta}_{k-1}$, and using a small learning rate during training. Alternative viable approaches include imposing trust-region constraints on the parameters of the form $\|\hat{\theta}_{k}-\hat{\theta}_{k-1}\|_{2}\leqslant\kappa$, or explicitly enforcing the trust-region constrain (4.1b). These latter options would be implemented through either a suitable Lagrangian relaxation to soft-penalties in the objective, or by drawing on recent results in constrained empirical risk minimization .

Enforcing the IGS constraint (4.1c) via an incremental Lyapunov function (cf. Proposition 3.3) is more challenging, as it must be enforced for all $x$ within a desired region of attraction. If such a Lyapunov function is known for the expert, then it can be used to only enforce stability constraints on trajectory data, an approximation/heuristic that is common in the constrained policy optimization literature (see for example ). However, if a Lyapunov certificate of IGS stability for the expert is not known, then options include (a) learning such a certificate for the expert from data, see for example , or (b) jointly optimizing over an IGS certificate and learned policy. Although this may be computationally challenging, alternating optimization schemes have been proposed and successfully applied in other contexts, see for example .

Fortunately, we note that empirically, explicit stability constraints seem not to be required. In the next subsection, we study the quantitative effects of enforcing the stability constraint (4.1c) for a linear system, for which a stability certificate is available, and for which the level of stability of the expert system can be quantitatively tuned. We observe that only when (a) the expert is nearly unstable, and (b) we are in a low-data regime, that a small difference in performance between stability-constrained and unconstrained algorithms occurs, suggesting that optimal policies are naturally stabilizing. Therefore, we simply omit constraint (4.1c) from our implementation and take care to ensure that sufficient data is provided to the IL algorithms to yield stabilizing policies.

2 Tuneable IGS System

In this experiment, we vary $p\in\{1,\dots,5\}$ to see the effect of $p$ on the final task goal error and imitation loss. We compare three different algorithms. BC is standard behavior cloning. CMILe is Algorithm 1 with the practical modifications as described above. DAgger is the imitation learning algorithm from Ross et al. . For each algorithm, we also consider a modification (indicated by the +IGS label) where policy imitation is augmented with a soft loss encoding the IGS constraint (4.1c). For all algorithms, we fix the number of trajectories $m$ from (5.1) to be $m=250$. The horizon length is $T=100$. The distribution $\mathcal{D}$ over initial condition is set as $N(0,I)$. We set the policy class $\Pi$ to be two layer MLPs with hidden width $64$ and $\tanh$ activations. Each algorithm minimizes the imitation loss using $300$ epochs of Adam with learning rate $0.01$ and batch size $512$. For all algorithms except BC, we use $E=25$ epochs with $\alpha=0.15$ (in DAgger’s notation, we set $\beta_{k}=0.85^{k}$), resulting in $10$ trajectories per epoch.

3 Unitree Laikago

We now study IL on the Unitree Laikago robot, an 18-dof quadruped with 3-dof of actuation per leg. We use PyBullet for our simulations. The goal of this experiment is to demonstrate, much like for the previous tuneable family of IGS systems, that increasing the stability of the underlying closed-loop expert decreases the sample-complexity of imitation learning. We do this qualitatively by studying a sideways walking task where the robot tracks a constant sideways linear velocity. By increasing the desired linear velocity, the resulting expert closed-loop becomes more unstable.

Our expert controller is a model-based predictive controller using a simplified center-of-mass dynamics as described in Di Carlo et al. . The stance and swing legs are controlled separately. The swing leg controller is based on a proportional-derivative (PD) controller. The stance leg controller solves for the desired contact forces to be applied at the foot using a finite-horizon constrained linear-quadratic optimal control problem; the linear model is computed from linearizing the center-of-mass dynamics. The desired contact forces are then converted to hip motor torques using the body Jacobian. More details about the expert controller can be found in the appendix.

We restrict our imitation learning to the stance leg controller, as it is significantly more complex than the swing leg controller. Furthermore, instead of randomizing over initial conditions, we inject randomization into the environment by subjecting the Laikago to a sequence of random push forces throughout the entire trajectory. We compare the performance of BC, CMILe, and CMILe+Agg; the CMILe+Agg algorithm is identical to CMILe, except that at epoch $k$, the data from previous epochs $j\in\{0,\dots,k-1\}$ is also used in training. DAgger is omitted for space reasons as its performance is comparable to CMILe+Agg.

We set the horizon length to $T=1000$, and featurized the robot state into a $14$-dimensional feature vector; the exact features are given in Appendix B.2. The output of the policy is a $12$-dimensional vector ($x,y,z$ contact forces for each of the $4$ legs). We used a policy class of two layer MLPs of hidden width $64$ with ReLU activations. For training, we ran $500$ epochs of Adam with a batch size of $512$ and step size of $0.001$. Furthermore, we tried to overcome the effect of overfitting in BC by using the following heuristic: we used $5\%$ of the training data as a holdout set, and we stopped training when either the holdout risk increased $h=50$ times or $500$ epochs were completed, whichever came first. To assess the effect of the number of samples on imitation learning, we vary the number of rollouts per epoch $S\in\{1,\dots,5\}$. For CMILe and CMILe+Agg, we fix $\alpha=0.3$ and $E=12$. We provide BC with $S\times E$ total trajectories.

Conclusions and Future Work

We showed that IGS-constrained IL algorithms allow for a granular connection between the stability properties of an underlying expert system and the resulting sample-complexity of an IL task. Our future work will focus on two complementary directions. First, CMILe and DAgger significantly outperform BC in our experiments, but our bounds are not yet able to capture this: future work will look to close this gap. Second, although our focus in this paper has been on imitation learning, we have developed a general framework for reasoning about learning over trajectories in continuous state and action spaces. We will look to apply our framework in other settings, such as safe exploration and model-based reinforcement learning.

Acknowledgements

The authors would like to thank Vikas Sindhwani, Sumeet Singh, Andy Zeng, and Lisa Zhao for their valuable comments and suggestions. NM is generously supported by NSF award CPS-2038873, NSF CAREER award ECCS-2045834, and a Google Research Scholar Award.

References

Appendix A Stability Study

for $\gamma\in(0,1)$ and $\varepsilon>0$. Note that by rewriting this equation as

Then note that we can rewrite the Lyapunov equation (A.1) as

for $Q=(1-\gamma^{2})X+\varepsilon I$. To that end, we suggest solving the Lyapunov equation

with $Q=(1-\gamma^{2})P_{\star}+\varepsilon I$, for a small $\varepsilon>0$ and $P_{\star}$ the solution to the DARE, so as to obtain a a Lyapunov certificate with similar convergence properties to that of the expert, where $\varepsilon$ trades off between how small $\gamma$ can be and how robust the resulting Lyapunov function is to mismatches between the learned policy and the expert policy. We note that the existence of solutions to this Lyapunov equation are guaranteed by continuity of the solution of the Lyapunov equation and that the solution to the DARE $P_{\star}$ is the maximizing solution among symmetric solutions .

A.2 Experimental Results

We study the effects of explicitly constraining the played policies $\{\pi_{i}\}_{i=1}^{E}$ to be IGS through the use of Lyapunov certificates. The use of Lyapunov constraints to enforce incremental stability was studied in Section G.1 of Boffi et al. , with the main takeaway being that systems satisfying suitable exponential Lyapunov conditions, in particular those certifying exponential input-to-state-stability, are also exponentially IGS (i.e., satisfy $a=b=1$) on a compact set of initial conditions and bounded inputs.

We drew initial conditions according to the distribution $N(0,4)$. We used a policy parameterized by a two-hidden-layer feed-forward neural network with ReLU activations. Each hidden layer in this network had a width of 64 neurons. For both BC and CMILe without stability constraints, we train the policy for 500 epochs; for CMILe with stability constraints, we trained policies for 1000 epochs. All neural networks were optimized with the Adam optimizer with a learning rate of $0.01$.

Appendix B Laikago Experimental Details

The expert controller contains multiple components: the swing controller, the stance controller, and the gait generator. The gait generator uses the clock source to generate a desired gait pattern, where a pair of diagonal legs are synchronized and are out of phase with the other pair. Throughout our experiments, we fixed the gait to be a trotting gait. The swing controller generates the aerial trajectories of the feet when they lift up and controls the landing positions based on the desired moving speed. The stance leg controller is based on model predictive control (MPC) using centroidal dynamics . Recall that in our experiments, we only perform imitation learning for the stance leg controller.

We now describe the centroidal dynamics model. We treat the whole robot as a single rigid body, and assume that the inertia contribution from the leg movements is negligible (cf. Figure 3). With these assumptions, the system dynamics can be simply written using the Newton-Euler equations:

where $\mathbf{x}=(x,y,z,\Phi,\Theta,\Psi)$ denotes the center of mass (CoM) translation and rotation, $\mathbf{F}_{i}=(f_{x},f_{y},f_{z})_{i}$ is the contact force applied on the $i$-th foot (set to zero if the $i$-th foot is not in contact with the ground), and $\mathbf{r}_{i}$ is the displacement from the CoM to the contact point. We used the same $Z-Y-X$ Euler angle conventions in to represent the CoM rotation. Since the robot operates in a regime where its base is close to flat, singularity from the Euler angle representation is not an concern.

The MPC module solves an optimization problem over a finite horizon $H$ to track a desired pose and velocity $\mathbf{q}=(\mathbf{x},\mathbf{\dot{x}})$. The system dynamics can be linearized around the desired state and discretized:

where $\mathbf{u}_{t}=(\mathbf{F}_{1,t},\mathbf{F}_{2,t},\mathbf{F}_{3,t},\mathbf{F}_{4,t})$, is the concatenated force vectors from all feet. We then formally write the optimization target:

where we used a diagonal $\mathbf{Q}$ and $\mathbf{R}$ matrix with the weights detailed in . At runtime, we apply the feet contact forces from the first step by converting them to motor torques using the Jacobian matrix.

B.2 Featurization

The inputs to the MPC algorithm is a 28 dimension vector $(\mathbf{q},\mathbf{q}^{d},\mathbf{c},\mathbf{r})$, where $\mathbf{c}$ is a binary vector indicating feet contact states, and $\mathbf{r}=(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3},\mathbf{r}_{4})$ represent the relative displacement from the CoM to each feet. This representation contains redundant information, since one can infer the body height $z$ from the contact state, and local feet displacements when the quadruped is walking on flat ground. Also, since in our experiments the desired pose and speed of the robot are fixed (moving along $y$ direction at constant speed without body rotation), they are not passed as inputs to the imitation policy. Furthermore, the current body linear velocities of the robot are omitted from the inputs as well, since they are not directly measurable on a legged robots without motion capture systems or state estimators. As a result, the inputs to the imitation policy are condensed to a 14 dimensional vector $(\Phi,\Theta,\mathbf{r\cdot c})$, i.e., the roll, pitch angle of the CoM, and the contact state masked feet positions.

Appendix C Incremental Gain Stability Proofs

We first prove a simple proposition which we use repeatedly.

We now derive some basic consequences of the definition of incremental gain stability. The following helper proposition will be useful for what follows.

For any $a\in[1,\infty)$ and integers $T_{1},T_{2}$ satisfying $0\leqslant T_{1}\leqslant T_{2}$, we have:

Let the index set $I\subseteq\{T_{1},\dots,T_{2}\}$ be defined as:

By Hölder’s inequality, since $a\in[1,\infty)$,

Next, we compare the autonomous trajectories between two different initial conditions (both trajectories are not driven by any input).

Now, suppose that $\lVert\Delta_{t}\rVert_{X}\leqslant 1$. By a similar argument:

Combining these inequalities yields the desired inequality (C.1):

Now we turn to (C.2). By Proposition C.2 and $(a,b,\Psi)$-incremental-gain-stability, we have:

The next result compares two trajectories starting from the same initial condition, but one being driven by an input sequence $\{u_{t}\}$ whereas the other is autonomous.

By Proposition C.2, the fact that $\Delta_{0}=0$, and $(a,b,\Psi)$-incremental-gain-stability,

Above, the last equality follows from Proposition C.1. The claimed inequality now follows. ∎

C.2 Proof of Proposition 3.3

Now define $V_{t}:=V(x_{t},y_{t})$. Then for $t\in\{0,\dots,T-1\}$, by the assumed inequality (3.6),

Next, by Proposition C.1, since $\alpha_{0}\in[1,a]$:

Appendix D Examples of Incremental Gain Stability Proofs

Recall the following definition of autonomously contracting in Proposition 3.4, which we duplicate below for convenience.

Consider the dynamics $x_{t+1}=f(x_{t},u_{t})$. We say that $f$ is autonomously contracting if there exists a positive definite metric $M(x)$ and a scalar $\gamma\in(0,1)$ such that:

For what follows, let $d_{M}$ denote the geodesic distance under the metric $M$:

The next result shows that the Euclidean norm lower and upper bounds the geodesic distance under $M$ as long as $M$ is uniformly bounded.

We now restate and prove Proposition 3.4. See 3.4

Fix initial conditions $\xi_{1},\xi_{2}$ and an input sequence $\{u_{t}\}_{t\geqslant 0}$. Consider two systems:

Above, (a) is triangle inequality, (b) follows from Proposition D.2 and Proposition D.3, and (c) follows from the Lipschitz assumption. Now unroll this recursion, to yield for all $t\geqslant 0$

Now dividing both sides by $\sqrt{\underline{\mu}}$ and summing the left hand side,

D.2 Scalar Example with p∈(0,∞)𝑝0p\in(0,\infty)

Recall we are interested in the family of systems:

Let us assume wlog that $y\geqslant|x|$, otherwise we can swap $x,y$ by considering $z(-y,-x)=z(x,y)$ instead. Now, we have

On the other hand, if $y\leqslant 1$, then

In this case, we have $z(x,y)=z(-y,-x)$. By swapping $x,y$, we reduce to Case 1 where we know that (D.1) already holds. ∎

Next, we show that the sign of $x-y$ is preserved under a perturbation $x-y-\eta(h(x)-h(y))$, as long as $\eta\geqslant 0$ is small enough.

Let $\eta$ satisfy $0\leqslant\eta<\frac{4}{5+p}$. We have that:

Observe that (D.2) holds trivially when $x=y$. Furthermore, when $y>0$:

Clearly $h^{\prime}(x)\geqslant 0$ when $x>0$. Furthermore, one can check that $\sup_{x\geqslant 0}\frac{x}{(1+x)^{2}}=1/4$. Therefore:

By assumption, we have that $1-\eta(1+(p+1)/4)>0$ and hence $z(x,y)<0$. This shows that $\operatorname*{sgn}(z(x,y))=\operatorname*{sgn}(x-y)$, and hence (D.2) holds in this case.

Hence, we have $\operatorname*{sgn}(z(x,y))=\operatorname*{sgn}(x-y)$, and hence (D.2) holds in this case.

Because $\operatorname*{sgn}(x)$ is an element of $\partial|x|$, by convexity of $|\cdot|$,

Above, (a) is Proposition D.5 and (b) is Proposition D.4. ∎

Note that Proposition 3.5 is an immediate consequence of Proposition D.6 with Proposition 3.3.

Appendix E Proof of Theorem 4.3 and Theorem 4.4

Our main tool will be the following uniform convergence result.

Next, define the following Rademacher complexity for the policy class $\Pi$:

Now fix a data generating policy $\pi_{d}\in\Pi(a,1,\Psi)$ and goal policy $\pi_{g}\in\Pi$. Furthermore, let $\xi_{1},\dots,\xi_{m}$ be drawn i.i.d. from $\mathcal{D}$. With probability at least $1-\delta$ (over $\xi_{1},\dots,\xi_{m}$), we have:

This follows from standard uniform convergence results, see e.g., Wainwright . ∎

Under Assumption 4.1 and Assumption 4.2, we have that:

Let $\pi_{d}\in\Pi(a,1,\Psi)$ and $\pi_{1},\pi_{2}\in\Pi$. Since $\Delta_{\pi_{1},\pi_{2}}(0)=0$ and $\Delta_{\pi_{1},\pi_{2}}$ is $L_{\Delta}$-Lipschitz:

Above, the last inequality follows from Proposition C.3. ∎

We now give a bound on the Rademacher complexity $\mathcal{R}_{m}(\Pi)$.

Fix an $x$ and $\theta_{1},\theta_{2}$. Since $\pi(0,\theta)=0$ for all $\theta$, by repeated applications of Taylor’s theorem:

Now supposing $\lVert x\rVert_{2}\leqslant\zeta B_{0}^{\alpha_{0}}$ and $\lVert\theta_{i}\rVert_{2}\leqslant B_{\theta}$ for $i\in\{1,2\}$, then:

The calculation above shows that for all $\pi_{1},\pi_{2}\in\Pi$,

Thus, for every $\varepsilon>0$, letting $L_{\Pi}:=2\zeta B_{0}^{\alpha_{0}}B_{\theta}L_{\partial^{2}\pi}T^{1-1/a}$, we have the following upper bound on the covering number:

Therefore by Dudley’s entropy integral (cf. Wainwright ):

The last inequality above follows from the numerical estimate:

Let $\operatorname*{conv}(\Pi)$ denote the convex hull of the policy class $\Pi$, and let

The following auxiliary proposition shows that the Rademacher complexity of the policy class $\operatorname*{conv}(\Pi)$ can be analyzed nearly identically to the Rademacher complexity of the original class $\Pi$.

Suppose that the policy class $\Pi$ is uniformly bounded, i.e.,

E.2 Proof

We first state our main meta-theorem, from which we deduce our rates.

Suppose that Assumption 4.1 and Assumption 4.2 hold. Suppose that $E\leqslant m$ divides $m$. Define $\Gamma(m,E,\delta)$ as:

Fix a $\delta\in(0,1)$. Assume that for all $k\in\{0,\dots,E-2\}$:

For $k\in\{0,\dots,E-1\}$, define $\beta_{k}(m,E,\delta)$ as:

With probability at least $1-\delta$ (over $\{\xi_{i}^{k}\}_{i=1,k=0}^{m/E,E-1}$ drawn i.i.d. from $\mathcal{D}$), we have that the following inequalities simultaneously hold for the policies $\pi_{1},\dots,\pi_{E}$ produced by Algorithm 1:

We first use induction on $k$ to show that:

Therefore, by our assumption that $\Delta_{\hat{\pi}_{0},\pi_{\star}}$ is $L_{\Delta}$-Lipschitz:

Now using the observation that $\Delta_{\pi_{1},\pi_{\star}}(x)=\alpha\Delta_{\hat{\pi}_{0},\pi_{\star}}(x)$, we write:

Therefore, on the event $\mathcal{E}_{0}$,

The first inequality above uses Jensen’s inequality to move the expectation inside $x\mapsto x^{1/a}$.

Furthermore we note that on $\mathcal{E}_{k}$, it holds that:

Above, (a) follows from (E.7) and (b) follows from (E.6).

Our remaining task is to show that on $\mathcal{E}_{0:k}$ we have:

We proceed with a similar argument as in the base case. We first write:

Therefore since by assumption $\Delta_{\pi_{k+1},\pi_{\star}}$ is $L_{\Delta}$-Lipschitz:

Combining this inequality with the inequality above,

Taking expectations of both sides, applying (E.6), and using Jensen’s inequality, we obtain

where the first inequality (a) follows from (E.7), and the second inequality (b) from $\hat{\pi}_{k}$ being feasible to the constrained optimization problem (4.1).

where (a) follows from (E.7), (b) from using $\pi_{k}$ as a feasible point for optimization problem (4.1) and optimality of $\hat{\pi}_{k}$, (c) from another application (E.7), and (d) follows from (E.6).

This finishes the inductive step. Thus we conclude that on the event $\mathcal{E}_{0:E-2}$, which occurs with probability at least $1-\frac{E-1}{E}\delta$, we have that for $k=1,\dots,E-1$,

We now assume $k=E-1$ and that the event $\mathcal{E}_{0:E-2}$ holds. On this event:

Furthermore we note that on $\mathcal{E}_{E-1}$, it holds that:

Therefore we can bound $c_{E-1}$ on $\mathcal{E}_{E-1}$ by:

Therefore since $\Delta_{\pi_{E},\pi_{\star}}$ is $L_{\Delta}$-Lipschitz by assumption,

Furthermore, it is straightforward to check that:

Combining this inequality with (E.13), taking expectations and applying Jensen’s inequality:

With Theorem E.5 in place, we now turn to the proof of our main results, which are immediate consequences of Theorem E.5. We first restate and prove Theorem 4.3. See 4.3

Theorem E.5 states that if $\Gamma(m,1,\delta)\leqslant 1$, then:

To complete the proof we simply need to bound $\Gamma(m,1,\delta)$, which has the form:

From Proposition E.2 and Proposition E.3, we have that:

We now restate and prove Theorem 4.4. See 4.4

We first bound $\Gamma(m,E,\delta)$, which has the form:

From Proposition E.2 and Proposition E.3, we have that:

Setting $\delta=e^{-q}$, this yields the bound:

We choose $m$ and $c_{k}$ such that (a) $\Gamma(m,E,\delta)\leqslant 1/2$ and (b) $c_{k}\leqslant\Gamma(m,E,\delta)$ for $k\in\{1,\dots,E-2\}$, which leads to the constraint (4.3) for (a) and the constraint (4.4) for (b).

In preparation to apply Theorem E.5, we use our assumptions to show:

Since $E\geqslant\frac{1}{\alpha}\log(\frac{1}{\alpha})$, then $\frac{(1-\alpha)^{E}}{\alpha}\leqslant\exp(-\alpha E)/\alpha\leqslant 1$. Furthermore, since we also assume that $\alpha\leqslant\min\left\{\frac{1}{2},\frac{1}{L_{\Delta}\gamma T^{1-1/a}}\right\}$, then

Now we proceed to bound $\beta_{E-1}(m,E,\delta)$:

Combining this bound with the inequalities for $(1-\alpha)^{E}$ yields (E.19).

We now apply Theorem E.5 with (E.19) and (E.20):

Since $E\Gamma(m,E,\delta)^{1/a}\leqslant 1$ by (4.3), we have: