Learning Stability Certificates from Data

Introduction

A fundamental barrier to widespread deployment of reinforcement learning policies on real robots is the lack of formal safety and stability guarantees. While much research has focused on how to train control policies for complex systems, considerably less emphasis has been placed on verifying stability for the resulting closed-loop system. Without any a-priori guarantees, practitioners will be hesitant to deploy learned solutions in the real world regardless of performance in simulation.

Many powerful tools have been developed in nonlinear control theory to address the safety and stability of systems with known dynamics. The most well-known technique is the construction of a Lyapunov function to demonstrate asymptotic stability of a system with respect to an equilibrium point. Similarly, barrier functions are used to show set-invariance, which has been widely used in safety-critical applications to prove that a system does not exit a desired safe set. Contraction analysis provides an alternative view of stability, applicable to many problems in nonlinear control and robotics, by considering the convergence of trajectories towards each other rather than to an equilibrium point. The unifying theme among these tools is the construction of a certificate function (the Lyapunov/barrier function or contraction metric) that proves a given desirable property for the system of interest. These certificates have strong converse results , which imply the existence of a certificate function if the desired property does hold, and can also be used for controller synthesis .

The main obstacle for producing certificate functions in modern robotics and reinforcement learning is that existing synthesis and verification tools such as sum-of-squares (SOS) optimization or SMT solvers typically assume the dynamics can be written down analytically in closed form. Furthermore, the functions of interest are often constrained to lie in restrictive classes such as polynomial basis functions of fixed degree. This presents a serious hurdle in modern robotics, where (a) sophisticated physics simulators are widely used to model complex environments and (b) control policies are often represented with complex deep neural networks. Finally, both SOS optimization and formal verification tools are computationally intensive, thus limiting their applicability.

To avoid these limitations, recent approaches have proposed to treat certificate synthesis as a machine learning problem, and train powerful function approximators such as deep neural networks and reproducing kernel Hilbert space (RKHS) predictors on trajectory data collected from a dynamical system . The general strategy is to enforce the desired certificate condition (e.g. the Lie derivative of a function $V$ should be negative) along collected samples. Empirically, this has been shown to be quite effective, and the learned certificate often generalizes well outside of the training data. However, a deeper theoretical understanding of when and why this approach works is missing.

Consider Figure 1, where a contraction metric, which certifies pairwise convergence of trajectories, is learned from rollouts of a damped Van der Pol oscillator. Regions of the state space for which the learned metric is not contracting are shown as a function of the number of trajectories $n$. While the size of the violating regions appears to shrink as $n$ increases, Figure 1 raises many questions. How much data does one need to collect so that the violating regions cover at most a prescribed fraction of the relevant state space? Is the learning consistent, i.e. do the regions vanish as $n\to\infty$?

Related Work

Prior research generally focuses on learning certificates for a fixed system from trajectories, or on using certificate conditions as regularizers when learning models for control.

Giesl et al. propose to learn a Lyapunov function from noisy trajectories using a specific reproducing kernel. Their algorithm first fits a dynamics model from data, and then uses interpolation to construct a Lyapunov function from the learned model. The authors prove $L_{\infty}$ convergence results on the Lie derivative of the constructed Lyapunov function compared to the ground truth, with rates depending on a dense cover of the state space.

Our work circumvents this two-step identification procedure by directly analyzing the generalization error of a Lyapunov function learned by enforcing derivative conditions along the training data. Many other authors have proposed similar approaches. Kenanian et al. show how to estimate the joint spectral radius of a switched linear system by learning a common quadratic Lyapunov function directly from data. Their analysis heavily exploits properties of linear systems. Chen et al. study how to learn a quadratic Lyapunov function for piecewise affine systems in feedback with a neural network controller. Richards et al. use a sum-of-squares neural network representation to learn the largest region of attraction of a nonlinear system. Manek and Kolter jointly train a neural network model and Lyapunov function. Neither Richards et al. nor Manek and Kolter provide formal guarantees that the learned Lyapunov function will generalize to new trajectories. Both Chang et al. and Ravanbakhsh and Sankaranarayanan propose to use ideas from formal verification to falsify the validity of a learned candidate Lyapunov function. A significant limitation is the requirement of access to the true dynamics.

In many of these works, the Lie derivative constraint that defines a Lyapunov function is relaxed to a soft constraint, so that first-order gradient methods can be used for optimization. We note that our generalization analysis can be modified to handle soft constraints in a straightforward manner.

Barrier functions are relaxations of Lyapunov functions that demonstrate invariance of a subset of the state space. Recently, many authors have proposed to use and learn barrier functions from data for safety-critical applications. Taylor et al. assume a control barrier function (CBF) is valid for both a nominal and unknown system model, and use the CBF to guide safe learning of the unknown system dynamics. More closely related to our work, Robey et al. learn a CBF for a known nonlinear dynamical system from expert demonstrations, and use Lipschitz arguments to extend the validity of the CBF beyond the training data. Jin et al. propose to jointly learn a Lyapunov, barrier, and a policy function from data. They also prove validity of the learned certificates using Lipschitz arguments.

The literature on learning contraction metrics from data is more sparse. In an imitation learning context, Sindhwani et al. propose to learn a vector field from demonstrations that satisfies contraction in the identity metric. The authors parameterize the vector field as a vector-valued reproducing kernel. Khadir et al. also learn a vector field from demonstrations by using sum-of-squares to enforce contraction. They argue by smoothness considerations that the learned vector field actually contracts in a tube around the demonstration trajectories. We note that in both these works, the metric is held fixed and is assumed to be known. Singh et al. jointly learn a model and a control contraction metric from data, and show empirically that using contraction as a regularizer in model learning can lead to better sample efficiency when learning to control. We leave studying the generalization properties of jointly learning an explicit model and a contraction metric to future work.

Our generalization bounds are similar in spirit to those provided for random convex programs (RCPs) . Random convex programming is concerned with approximating solutions to convex programs with an infinite number of constraints. Such infinitely-constrained problems are approximated by drawing $n$ i.i.d. samples from a distribution $\nu$ over the constraint parameters and enforcing constraints on samples. One can then show that the probability that a new sample from $\nu$ violates the constraint for the approximate solution scales as $O(d/n)$ where $d$ is the number of decision variables. Our results can be viewed as generalizing these bounds beyond convex programs, though our constants are less sharp. In our experiments, we use the RCP bound for numerically computing generalization bounds when the problem is convex.

Learning Certificates Framework

As we describe below, through suitable choices of the function $h$, equation (3.1) can be used to enforce various defining conditions for certificates such as Lyapunov functions and contraction metrics. We note that our framework can be modified to allow for more derivatives of $V$, including higher order derivatives and also time derivatives for handling time-varying dynamics.

We study the following optimization problem for searching for a solution to (3.1):

Here, $\gamma>0$ is a positive margin value which will allow us to generalize the behavior of $V$ on $h$ outside of the sampled data. In practice, we often solve (3.2) with a cost term on $V$ such as its norm. Let $\hat{V}_{n}\in\mathcal{V}$ denote a solution to (3.2), assuming one exists. We quantify the generalization of $\hat{V}_{n}$ by the probability of violation over trajectories starting from $\xi\sim\mathcal{D}$:

1.2 Contraction metrics

A system is said to be contracting in a region $D$ with rate $\alpha$ if there exists a uniformly positive definite Riemannian metric $M(x)$ such that $\frac{\partial f}{\partial x}(x)^{\mathsf{T}}M(x)+M(x)\frac{\partial f}{\partial x}(x)+\dot{M}(x)\preceq-2\alpha M(x)$ for $x\in D$ . Given knowledge of $\frac{\partial f}{\partial x}$, this condition fits into our framework by taking $h\left(x,\dot{x},M(x),\frac{\partial f}{\partial x}(x)\right)=\lambda_{\max}\left(\frac{\partial f}{\partial x}(x)^{\mathsf{T}}M(x)+M(x)\frac{\partial f}{\partial x}(x)+\dot{M}(x)+2\alpha M(x)\right)$.

We can enforce (3.4) by directly sampling trajectories on $\mathcal{T}D$, by exploiting that the variational dynamics obeyed by $\delta x(t)$ is identical to the local linearization of $f$ around $x(t)$. Specifically, we sample pairs of initial conditions $x^{(1)}_{0}$ and $x^{(2)}_{0}=x^{(1)}_{0}+\delta x_{0}$ for some small perturbation $\delta x_{0}$. Numerical differentiation of $x^{(1)}(t)$ and $\delta x(t)=x^{(1)}(t)-x^{(2)}(t)$ provides access to $\dot{x}^{(1)}=f(x^{(1)})$ and $\delta\dot{x}(t)=\frac{\partial f}{\partial x}\delta x(t)$, which then allows us to evaluate (3.4) along system trajectories.

Generalization Error Results

We first define the notion of stability we will assume. Recall that $X$ is the set containing sample initial conditions, and $T$ is the interval over which our trajectories evolve.

Note that contraction implies Assumption 4.1, so that contracting systems are also covered in this setting. Next, we make some regularity assumptions on the function class $\mathcal{V}$.

We assume there exist finite constants $B_{V}$, $B_{\nabla V}$ such that $\sup_{V\in\mathcal{V}}\sup_{x\in S}\lVert V(x)\rVert\leq B_{V}$ and $\sup_{V\in\mathcal{V}}\sup_{x\in S}\left\lVert\frac{\partial V}{\partial x}(x)\right\rVert\leq B_{\nabla V}$.

Given Assumptions 4.1–4.2, we define $B_{h}$ (resp. $L_{h}$) to be an upper bound on $|h(x,f(x),V,\frac{\partial V}{\partial x})|$ (resp. the Lipschitz constant of $(V,\frac{\partial V}{\partial X})\mapsto h(x,f(x),V,\frac{\partial V}{\partial X})$) over $x\in S$, $|V|\leq B_{V}$ and $\lVert\frac{\partial V}{\partial x}\rVert\leq B_{\nabla V}$. Note that both $B_{h}$ and $L_{h}$ are guaranteed to be finite by our assumptions.

Fix a $\delta\in(0,1)$. Assume that Assumption 4.1 and Assumption 4.2 hold. Suppose that the optimization problem (3.2) is feasible and let $\hat{V}_{n}$ denote a solution. The following statement holds with probability at least $1-\delta$ over the randomness of $\xi_{1},...,\xi_{n}$ drawn i.i.d. from $\mathcal{D}$:

We consider the following parametric representation:

Under Assumption 4.1, if problem (3.2) over the parametric function class (4.1) is feasible, then any solution $\hat{V}_{n}$ satisfies with probability at least $1-\delta$ over $\xi_{1},...,\xi_{n}$:

Here, $L_{g}:=\sup_{x\in S,\lVert\theta\rVert\leq B_{\theta}}\lVert\nabla_{\theta}g(x,\theta)\rVert$ and $L_{\nabla g}:=\sup_{x\in S,\lVert\theta\rVert\leq B_{\theta}}\lVert\frac{\partial^{2}g}{\partial\theta\partial x}(x,\theta)\rVert$.

Often times (4.1) is more structured. For instance, in sum-of-squares (SOS) optimization, we have:

Under Assumption 4.1, if problem (3.2) over the parametric linear function class (4.3) is feasible, then any solution $\hat{V}_{n}$ satisfies with probability at least $1-\delta$ over $\xi_{1},...,\xi_{n}$:

Here, $B_{m}:=\sup_{x\in S}\lVert m(x)\rVert$ and $B_{Dm}:=\sup_{x\in S}\lVert\frac{\partial m}{\partial x}(x)\rVert$.

In general, using prior knowledge about the system to add more structure and reduce the number of parameters of the certificate function (e.g. diagonal contraction metrics for positive systems) yields better generalization bounds.

2 Reproducing kernel Hilbert space function classes

We now consider the following non-parametric function class:

Here, $\phi(\cdot;\theta)$ is a nonlinear function and $\nu$ is a probability distribution over $\Theta$. This function class is a subset of the reproducing kernel Hilbert space (RKHS) defined by the kernel $k(x,y)=\int_{\Theta}\phi(x;\theta)\phi(y;\theta)\nu(\theta)\>d\theta$, and is dense in the RKHS as $B_{\alpha}\to\infty$ . We further assume that $\phi(x;\theta)$ is of the form $\phi(x;\theta)=\phi(\langle x,w\rangle+b)$ with $\phi$ differentiable and $\theta=(w,b)$. RKHSs of this type often arise naturally. For instance, Bochner’s theorem states that every translation invariant kernel can be expressed in this form.

Suppose that $|\phi|\leq 1$, $\phi$ is $L_{\phi}$-Lipschitz, $\phi$ is differentiable, $\phi^{\prime}$ is $L_{\phi^{\prime}}$-Lipschitz, and that $B_{\theta}:=\sup_{\theta\in\Theta}\lVert\theta\rVert$ is finite. Under Assumption 4.1, if problem (3.2) over the non-parametric class (4.5) is feasible, then any solution $\hat{V}_{n}$ satisfies with probability at least $1-\delta$ over $\xi_{1},...,\xi_{n}$:

where $\kappa:=\frac{B_{\alpha}^{2}L_{h}^{2}}{\gamma^{2}}\left((1+B_{\theta}L_{\phi})^{2}\log{n}+B_{\theta}^{2}(B_{S}+1)^{2}(L_{\phi}+B_{\theta}L_{\phi^{\prime}})^{2}p\right)$.

Global Stability Results

In this section, we show how the bounds from Section 4 can be translated into global results for the learned certificate functions. To facilitate our analysis, we assume the dynamics is incrementally stable. Incremental stability is implied by contraction, but is stronger than Lyapunov stability. Before stating the assumption, we say that $\beta(s,t)$ is a class $\mathcal{KL}$ function if for every $t$ the map $s\mapsto\beta(s,t)$ is a class $\mathcal{K}$ function and for every $s$ the map $t\mapsto\beta(s,t)$ is continuous and non-increasing.

There exists a class $\mathcal{KL}$ function $\beta$ such that for all $\xi_{1},\xi_{2}\in X$, $\lVert\varphi_{t}(\xi_{1})-\varphi_{t}(\xi_{2})\rVert\leq\beta(\lVert\xi_{1}-\xi_{2}\rVert,t)$ for all $t\in T$.

Equation (5.3) yields bounds of the form $V(\varphi_{t}(\xi))\leq V(\xi)e^{-\lambda t}+O(r(\varepsilon)h_{\beta}(t))$, where $h_{\beta}$ depends on the specific form of $\beta$. For example if $\beta(s,t)\leq Me^{-\alpha t}s$ for some $\alpha>\lambda$, then $h_{\beta}(t)=e^{-\lambda t}$. On the other hand, if we have the slower rate $\beta(s,t)\leq Ms/(t+1)$, then $h_{\beta}(t)=t^{2}e^{-\lambda t}$.

We note that Theorem 5.1 is conceptually similar to the results from Liu et al. , but incremental stability assumption dramatically simplifies the proof and enables us to make the constants explicit.

Theorem 5.2 is illustrated in Figure 1, which shows the structure of the violation set. Further details and exploration of the effect of $\eta$ can be found in Section A.5. In Section D we prove a result similar to Theorem 5.2 for metric learning with known dynamics.

Learning Certificates in Practice

We empirically study the generalization behavior of both learning Lyapunov functions and contraction metrics from trajectory data. We consider Lyapunov functions parameterized by $V(x)=x^{\mathsf{T}}(L(x)L(x)^{\mathsf{T}}+I)x$, where $L(x)$ is the (reshaped) value of a fully connected neural network with $\tanh$ activations of size $p\times h\times h\times p\cdot(2p)$, where $p$ is the state-dimension of $x$ and $h$ is the hidden width. For metric learning, we study a convex formulation via SOS programming. Each matrix element $M_{ij}(x)=\langle w_{ij},\phi(x)\rangle$ is given by a polynomial where $w_{ij}$ are the learned weights and $\phi(x)$ is a feature map of monomials in the state vector. In our experiments, we numerically estimate the generalization error of a learned certificate using a test set. We compute an upper confidence bound (UCB) of the estimate using the Chernoff inequality with $\delta=0.01$, as described by Langford . More experimental details are given in the appendix.

We learn a discrete-time Lyapunov function for a quadruped robot as it recovers from external forcing. We apply a random impulse force in the $(x,y)$ plane at time $t=0$ to the Minitaur quadruped environment in PyBullet , and use a hand-tuned PD controller to return the minitaur to a standing position. We train a discrete-time Lyapunov function in order to handle the discontinuities in the trajectories introduced by contact forces.

Figure 2 (LL) shows the result of this experiment. For the Lyapunov curve, the resulting model trained on $n$ trajectories is then validated using a $10000$ trajectory test set. The generalization error is the ratio of trajectories which violate the desired decrease condition for any step $k\in\{1,...,199\}$. We run $30$ trials and plot the 10/50/90-th percentile of the generalization UCB. With $n=50000$, the median generalization UCB is $1.11\%$. Since in practice a separate test set may not be available, we also compare to splitting the available training data into an actual training set of size $0.9n$ and a validation set of size $0.1n$. The model is trained on the actual training set, and a generalization UCB is calculated from the validation set. We run $30$ trials of this setup and plot the 10/50/90-th percentile in the Holdout curve. After $n=50000$, the median generalization UCB is $1.99\%$.

Conclusion

Our work shows that certificate functions can be efficiently learned from data, and raises many interesting questions for future work. Extending the results to handle both noisy state observations and process noise in the dynamics would allow for learning certificates in uncertain environments. Another interesting question is to establish bounds for joint learning of both the unknown dynamics and a certificate, which has shown to be effective in practice . Finally, lower bounds on the learning certificate problem would highlight the amount of conservatism introduced in our results.

Acknowledgements

The authors would like to thank Amir Ali Ahmadi, Brett Lopez, Alexander Robey, and Sumeet Singh for providing helpful feedback.

References

Appendix A Experiment Details

Here we give pseudocode for the metric learning algorithm used in the main text. Algorithm 1 is written for a parameterization $M_{w}(x)$ that yields a convex optimization problem, and where uniform positive definiteness may be enforced globally, such as via SOS matrix constraints. It may be readily relaxed to nonconvex parameterizations such as neural networks by using soft constraints and minimizing the loss using a variant of stochastic gradient descent. Uniform positive definiteness can be imposed along trajectories rather than globally.

A.2 Pendulum

We repeat this experiment for $30$ trials. For each trial we use a test set of size $1000$ to compute a UCB on the generalization error. The 10/50/90-th percentile of the UCBs are $0.459\%$, $0.459\%$, and $1.163\%$.

We also uniformly grid the set $\times$ with $40000$ points and check numerically how often the condition $\langle\nabla V_{\theta}(x),\dot{x}\rangle\leq-\gamma V_{\theta}(x)$ is violated. The 10/50/90-th percentile of these violations over $30$ trials are $0.168\%$, $1.548\%$, and $3.696\%$. We note that these numbers are higher than the generalization error because the set $\times$ contains points that are outside the flow starting from $\times$.

For our adaptive control experiments, the dynamics combined with the added disturbance are

with $f$ continuously differentiable and $Y(x,t)$ locally bounded in $x$ uniformly in $t$. Let $V$ be a twice continuously differentiable positive definite function such that $\langle\nabla V(x),f(x)\rangle\leq-\rho(x)$ for all $x$ for some continuously differentiable positive definite function $\rho(\cdot)$. Let $M$ be a positive definite matrix. Let $\hat{a}(t)$ be defined by the differential equation

in feedback with (A.2) drives $x\to 0$ and $\dot{x}\to 0$.

Let $\langle x,y\rangle_{M}=x^{\mathsf{T}}My$ and $\lVert x\rVert_{M}^{2}=\langle x,x\rangle_{M}$. Define the new candidate Lyapunov function

A.3 Minitaur

We collect $50000$ random training trajectories and $10000$ random test trajectories using the same distribution over the impulse force. For each trajectory, we step the simulator $200$ times at $dt=0.002$. The state dimension excluding the base orientation is $16$.

We use a hidden width of $h=40$ and minimize the loss $L_{n}(\theta)=\sum_{i=1}^{n}\sum_{k=1}^{199}\mathsf{ReLU}(V_{\theta}(e_{i}(k+1))-\rho V_{\theta}(e_{i}(k))-\gamma)+\lambda\lVert\theta\rVert^{2}$, setting $\lambda=0.01$. The loss is minimized for $1000$ epochs with Adam using a step size $10^{-3}$ with cosine decaySee https://www.tensorflow.org/api_docs/python/tf/compat/v1/train/cosine_decay. and a batch size of $1000$.

A.4 6-dimensional gradient system

We parametrize the metric via monomials up to degree two in the state variables. We enforce global positive definiteness $M(x)\geq\mu I$ via SOS matrix constraints and set $\mu=1$. We use a tolerance of $5\times 10^{-3}$ for the size of each perturbation $\delta x$. A pair of trajectories $\varphi_{t}(x_{1})$, $\varphi_{t}(x_{2})$ with $x_{2}=x_{1}+\delta x$ is considered to generate a trajectory $\delta x(t)=\varphi_{t}(x_{2})-\varphi_{t}(x_{1})$ if $\|\delta x(t)\|<10^{-2}$ for all $t$, so that a small overshoot is permitted. Pairs of trajectories not satisfying this requirement are discarded until the desired number of training samples is reached. The size of this overshoot parameter sets the maximum allowed $L$ where $M(x)\leq L$ for any metric learned, as we impose $M(x)\geq I$. In general, the overshoot with respect to the Euclidean norm is given by $\sqrt{\frac{L}{l}}$ for $lI\leq M(x)\leq LI$. In practice, we require the maximum bound on $\|\delta x(t)\|$ to be sufficiently small that the dynamics of $\varphi_{t}(x_{1})-\varphi_{t}(x_{2})$ well-approximates the variational system on the trajectory $\varphi_{t}(x_{1})$. To search for metrics with larger values of $\frac{L}{l}$, we can vary $\mu$, $\|\delta x\|$, and the maximum allowed $\|\delta x(t)\|$ while ensuring $\|\delta x(t)\|$ remains small throughout its entire trajectory.

Each trajectory is simulated until $T=2$ seconds with a timestep $dt=5\times 10^{-3}$. Because we are interested in convergence of the variational dynamics, we use a small time horizon. This generates a dataset of size $n\times 400\times 12$ where $n$ is the number of trajectories and $12$ is the dimension of the tangent bundle. We subsequently downsample and impose $25$ differential Lyapunov constraints along each trajectory. We search for a metric with a rate $\lambda=4$. Initial conditions are drawn uniformly from the ball of radius $r=3$.

The time derivatives $\dot{x}$ and $\delta\dot{x}$ are computed numerically by fitting a cubic spline to the corresponding trajectories and analytically differentiating the spline. $M(x)$ is found by minimizing $\|w\|^{2}$ where $w$ is a vector containing all parameters. The test set is of size $1000$ and each data point in Figure 2 (LR) was computed by averaging over $25$ independent draws of the training set. The RCP bound is obtained with a confidence of $\delta=.01$.

Consider the following optimization problem:

Fix any $\varepsilon\in$. Define $\beta(\varepsilon)$ as:

With probability at least $1-\beta(\varepsilon)$ over $\theta_{1},...,\theta_{n}$, we have that:

We use Theorem A.2 as follows. We fix a failure probability $\delta\in(0,1)$. We then numerically solve for $\varepsilon$ such that $\beta(\varepsilon)=\delta$.

To understand the scaling of $\varepsilon$ as a function of $n$ and $\delta$, despite the lack of closed form expression, consider the following. If $n\geq d$ then by a Chernoff bound on the CDF of a Binomial random variable (c.f. Section 5 of Calafiore ), we can derive the upper bound

A.5 Van der Pol

The study of the Van der Pol (VDP) oscillator was foundational to the development of nonlinear dynamics , and its global contraction properties have been analyzed algoithmically via SOS programming . We study the damped VDP to visualize the violation set for the metric condition categorized in Theorem 5.2. The dynamics of the damped Van der Pol are given by

We parameterize the metric via monomials up to degree four in the state variables. Similar to the gradient system, we set $M(x)\geq I$ and use a tolerance $\|\delta x\|\leq 5\times 10^{-3}$, $\|\delta x(t)\|\leq 10^{-2}$ for all $t\in[0,T]$. We use a timestep $dt=5\times 10^{-3}$ and simulate until a final time $T=3$ seconds. This generates a dateset of size $n\times 600\times 4$ where $4$ is the dimension of the tangent bundle and $n$ is the number of training trajectories.

We subsequently downsample and impose $50$ differential Lyapunov constraints along each trajectory with a rate $\lambda=3/4$. Initial conditions are drawn uniformly from a ball of radius $r=2$. The same techniques are used for numerical differentiation as for the gradient system.

Theorem 5.2 predicts that the size of the violation set will decrease as the rate tested for the metric condition $\eta\lambda$ with $0<\eta<1$ decreases, or as the number of training samples $n$ increases. In both Figure 1 and Figure 3, we use a uniform grid with $9\times 10^{4}$ points over $\times$ to test the metric condition $\frac{\partial f}{\partial x}(x)^{\mathsf{T}}M(x)+M(x)\frac{\partial f}{\partial x}(x)+\dot{M}(x)\leq-2\eta\lambda$ for the learned metric and the true dynamics.

In Figure 1, we plot the violation set for fixed $\eta=1$ as a function of the number of training samples. As discussed in the main text, the size of the violation set decreases, and it is pushed to the boundary of the sampled region as $n$ increases.

In Figure 3, we plot the violation set as a function of $\eta$. As $\eta$ decreases to zero, the size of the violation set decreases and is pushed to the boundary of the sampled region.

Appendix B Proofs for Section 4

Recall that Dudley’s inequality gives us the following estimate on $\mathcal{R}_{n}(\mathcal{V})$:

B.2 Proof of Theorem 4.3

A simple calculation shows that $\lVert V_{Q_{1}}-V_{Q_{2}}\rVert_{\mathcal{V}}\leq(B_{m}^{2}+2B_{Dm}B_{m})\lVert Q_{1}-Q_{2}\rVert$. Hence by (B.1),

Here, $c$ is an absolute constant. By integrating this estimate we arrive at the bound:

B.3 Proof of Theorem 4.4

Define the function classes $\mathcal{F}(B)$ and $\hat{\mathcal{F}}(B,\{\theta_{k}\})$ as:

The following is a simple modification of Theorem 3.2 from Rahimi and Recht , which also accounts for uniform approximation of the derivatives.

Furthermore, now assume that $\phi$ is differentiable and $\phi^{\prime}$ is $L_{\phi^{\prime}}$-Lipschitz. Then every $f\in\mathcal{F}(B)$ is differentiable with

Finally, with probability at least $1-2\delta$, there exists a $\hat{f}\in\hat{\mathcal{F}}(\lVert f\rVert_{\nu},\{\theta_{k}\}_{k=1}^{K})$ such that (B.5) holds and also:

Following the proof of Theorem 3.2 of Rahimi and Recht , we set $c_{k}=\frac{\alpha(\theta_{k})}{K\nu(\theta_{k})}$, and we define $v(\theta_{1},...,\theta_{K})=\lVert\hat{f}-f\rVert_{\infty}$. It is shown in that for all $\theta_{1},...,\theta_{K},\theta^{\prime}_{k}\in\Theta$:

Next, we control the expected value of $\nu$:

Above, the first inequality is a standard symmetrization argument where the Rademacher variables $\{\varepsilon_{k}\}$ are introduced, and the second inequality is the triangle inequality. We first bound $T_{1}$. Set $\psi_{k}(z):=c_{k}(\phi(z)-\phi(0))$. Clearly $\psi_{k}(0)=0$, and also $\psi_{k}$ is $c_{k}L_{\phi}$-Lipschitz. We bound $c_{k}L_{\phi}\leq\frac{\lVert f\rVert_{\nu}L}{K}$. Therefore by the contraction inequality for Rademacher complexities (Theorem 4.12 of Ledoux and Talagrand ) followed by Jensen’s inequality:

Furthermore we can bound $T_{2}\leq\frac{2\lVert f\rVert_{\nu}}{\sqrt{K}}$ by similar arguments. The claim (B.5) now follows by invoking McDiarmid’s inequality.

We note that (B.6) follows from a basic application of the dominated convergence theorem, since we have that:

Finally, we focus on the derivative condition (B.7). Let $g(\theta_{1},...,\theta_{K}):=\sup_{x\in X}\lVert\nabla\hat{f}(x)-\nabla f(x)\rVert$. By symmetrization we have:

We can now apply Theorem 3 of Maurer to conclude that:

Next, we have for all $\theta_{1},...,\theta_{K},\theta^{\prime}_{k}\in\Theta$:

The uniform bound on the derivatives (B.7) now follows from another application of McDiarmid’s inequality. ∎

We now turn to the proof of Theorem 4.4. Under the hypothesis of Lemma B.1, we have that for every $f\in\mathcal{F}(B)$:

We also have for any $\{\theta_{k}\}_{k=1}^{K}\subseteq\Theta$ and any $\hat{f}(x)=\sum_{k=1}^{K}c_{k}\phi(x;\theta_{k})$ with $\lVert c\rVert_{1}\leq B$,

Hence for any $\{\theta_{k}\}_{k=1}^{K}\subseteq\Theta$, the function class $\hat{\mathcal{F}}(B_{\alpha},\{\theta_{k}\}_{k=1}^{K})$ satisfies Assumption 4.2 with $B_{V}=B_{\alpha}$ and $B_{\nabla V}\leq B_{\alpha}B_{\theta}L_{\phi}$.

Let $f_{n}\in\mathcal{F}(B_{\alpha})$ denote a feasible solution to (3.2).

At this point, it may be tempting to use the probabilistic method in conjunction with Lemma B.1 to conclude that there exists a set of weights $\{\overline{\theta}_{k}\}_{k=1}^{K}$ such that there exists a $\hat{f}_{n}\in\hat{\mathcal{F}}(B_{\alpha},\{\overline{\theta}_{k}\}_{k=1}^{K})$ such that $(\hat{f}_{n},\nabla\hat{f}_{n})$ closely approximates $(f_{n},\nabla f_{n})$. This will not work however, since the function class $\hat{\mathcal{F}}(B_{\alpha},\{\overline{\theta}_{k}\}_{k=1}^{K})$ then becomes a function of the training data $\xi_{1},...,\xi_{n}$, and hence we would not be able to apply Lemma 4.1 to it.

To work around this, we need to draw the weights independently of $\xi_{1},...,\xi_{n}$. In particular, we set $K$ such that

where $c$ is an absolute constant and let $\{\theta^{*}_{k}\}_{k=1}^{K}$ be drawn i.i.d. from $\nu$.

By the definition of $L_{h}$, these two inequalities imply that

The result is that with probability at least $1-\delta$ over $\xi_{1},...,\xi_{n}$:

Letting $\hat{f}=\sum_{k=1}^{K}c_{k}\phi(x;\theta^{*}_{k})$ and $\hat{g}=\sum_{k=1}^{K}d_{k}\phi(x;\theta^{*}_{k})$, we have that $\lVert\hat{f}-\hat{g}\rVert_{\mathcal{V}}\leq(1+B_{\theta}L_{\phi})\lVert c-d\rVert_{1}$. Hence by (B.1),

Appendix C Proofs for Section 5

Now if $\mu(U)\leq\varepsilon$, then $\mu(U)=\frac{\mu_{\mathsf{Leb}}(U)}{\mu_{\mathsf{Leb}}(X)}\leq\varepsilon$ and hence $\mu_{\mathsf{Leb}}(U)\leq\varepsilon\mu_{\mathsf{Leb}}(X)$.

where $I(x;a,b)$ is the regularized incomplete beta function.

The solution to an LTV system is given by $x(t)=\Phi(t)x(0)$, where $\Phi(t)=\exp(\int_{0}^{t}A(\tau)\>d\tau)$ and $\Phi(t)$ is invertible for all $t$. This means there exists a non-zero $\xi$ such that $z=\Phi(t)\xi=\Phi(t)\frac{\xi}{\lVert\xi\rVert}\lVert\xi\rVert$. The claim now follows by taking $\alpha=1/\lVert\xi\rVert$. ∎

We therefore have the following chain of inequalities:

C.2 Proof of Theorem 5.2

We begin with a simple lemma, which shows that if a system evolving on Euclidean space is contracting in the metric $M(x,t)$, then the corresponding prolongated system on the tangent bundle will be contracting in a block-diagonal metric.

Consider the differential dynamics $\delta\dot{x}=\frac{\partial f}{\partial x}(x,t)\delta x$. This system has Jacobian

Consideration of the differential Lyapunov function

shows that $V$ decreases exponentially by contraction of $f(x,t)$ in the metric $M(x,t)$ and hence that the virtual dynamics are contracting. Let $\Theta(x,t)$ be such that $\Theta^{\mathsf{T}}\Theta=M$. The metric transformation

for $\epsilon>0$ leads to the generalized Jacobian

which shows that the prolongated system is contracting over any compact domain for $\epsilon$ sufficiently small. In particular, for contraction with rate $\eta\gamma$ for $0<\eta<1$, we may set

Furthermore, note that the metric transformation $\Theta^{\prime}$ corresponds to the block-diagonal metric

Now with this expression in hand, we choose $r_{b}$ such that

By $2$-homogeneity of the inequality above, we may divide by $\alpha^{2}$ on both sides to conclude:

Appendix D Known dynamics

In Section 5, we assume access only to trajectories. Here we prove a simple proposition under the assumption that the dynamics is known, so that the defining metric condition for contraction can be sampled directly.

with an identical bound for the transpose. Now let $\langle\nabla M,\dot{x}\rangle$ denote the tensor contraction $\langle\nabla M,\dot{x}\rangle_{ij}=\langle\nabla M_{ij},\dot{x}\rangle$. Then,

Furthermore, $M_{g}-M_{b}\leq L_{M}\|\xi_{b}-\xi_{g}\|$. Putting these bounds together, we find that