Single Trajectory Nonparametric Learning of Nonlinear Dynamics

Introduction

Estimation of models (1) and (2) remains relatively poorly understood when the data is not i.i.d., with existing results being limited to when the function $f_{\star}$ is known to belong to certain parametric classes. In terms of parameter recovery, the LSE converges at a rate of $T^{-1/2}$ for stable linear autoregressive systems $f_{\star}(x_{t})=A_{\star}x_{t}$ (Simchowitz et al., 2018; Sarkar and Rakhlin, 2019; Jedra and Proutiere, 2020). The same rate can also be achieved for linear systems with more general input-output behavior (Oymak and Ozay, 2019; Tsiamis and Pappas, 2019). Moving to nonlinear models, recursive and gradient type algorithms can be shown to converge at a rate of $T^{-1/2}$ for the generalized linear model $f_{\star}(x_{t})=\phi(A_{\star}x_{t})$, where $\phi$ is a known Lipschitz link function (Foster et al., 2020; Sattar and Oymak, 2020; Jain et al., 2021). In this paper, we significantly generalize these results and provide rate-optimal error bounds for nonparametric function classes in terms of their metric entropies. Our approach leverages recently developed information-theoretic tools (Russo and Zou, 2019; Xu and Raginsky, 2017) and the notion of offset complexity (Rakhlin and Sridharan, 2014; Liang et al., 2015).

Given this, the learning task is to produce an estimate $\hat{f}$ of the model $f_{\star}$, which is evaluated in terms of the expected Euclidean $2$-norm error:

The expectation (3) is computed with respect to the randomness in the algorithm and the random variable $\xi\sim\nu$ which is independent of all other randomness; here $\nu$ is the uniform mixture over $(x_{0},\dots,x_{T-1})$. That is, $\xi$ has the same distribution as random variable $x_{\tau}$, where the index $\tau$ is drawn uniformly at random over $\{0,\dots,T-1\}$ and is independent of $(x_{0},\dots,x_{T-1})$. If the process (1) is stationary with invariant measure $\nu$, this reduces to the assumption that $\xi$ is drawn from the invariant measure of the process.

In the sequel, we analyze the nonparametric least squares estimator (LSE) $\hat{f}$ of $f_{\star}$, given below. Namely, we assume that the learner can compute

1 Contributions

2 Further Related Work

Estimation of models of the form (1) and (2) has a rich history in statistics and system identification (Ljung, 1999). Preceding the recent body of work mentioned in the introduction, asymptotically optimal rates for linear stochastic models have been available for some time (Mann and Wald, 1943; Lai and Wei, 1982). Similarly, there is a well-established theory of rate-optimal identification for nonlinear parametric models under various identifiability-type conditions, both in the i.i.d. setting (Van der Vaart, 2000) and under more general assumptions (Le Cam, 2012).

Perhaps the main motivator for the recent line of work emphasizing nonasymptotic estimation bounds is that these bounds are applicable downstream in control and reinforcement learning pipelines. Regression estimates for linear stochastic systems have been key to understanding both online and offline reinforcement learning in the linear quadratic regulator (Dean et al., 2020; Mania et al., 2019) and can be shown to lead to optimal regret rates (Simchowitz and Foster, 2020; Ziemann and Sandberg, 2022). Extending our understanding of the interaction between learning and control beyond linear-in-the-parameters models (Kakade et al., 2020; Boffi et al., 2021; Lale et al., 2021) inevitably requires new analyses of learning in dynamical systems. This also motivates the present work in that we provide nonasymptotic and counterfactual control of the LSE’s estimation error for more general nonlinear and nonparametric models.

Another related field is that of general statistical learning for dependent data, see Agarwal and Duchi (2012), Kuznetsov and Mohri (2017) and the references therein. These works provide generalization bounds for general loss functions and $\beta$-mixing processes $(x_{t},y_{t})$. The assumption of $\beta$-mixing processes has also previously been exploited in parametric identification by Vidyasagar and Karandikar (2006). By contrast, our emphasis on regression over general learning is motivated by downstream applications in learning-enabled control, where one first learns a model used to design a controller. Necessarily then, this work builds on a rich line of work in nonparametric regression for the i.i.d. setting, see chapters 13 and 14 of Wainwright (2019) and the references therein. Although there has been some work on the dependent setting, the error in this line of work is typically computed with respect to the design points which is not suitable for the counterfactual reasoning that is key to control (Baraud et al., 2001).

We also draw inspiration from the recent line of work on information-theoretic generalization bounds (Xu and Raginsky, 2017; Russo and Zou, 2019; Bu et al., 2020). There are interesting refinements and variations of this theory using for instance conditional mutual information (Steinke and Zakynthinou, 2020), or Wasserstein distance (Gálvez et al., 2021). However, these more recent bounds rely more explicitly on the tenzorization properties of information measures under i.i.d. data than the earlier work of Russo and Zou (2019) and Xu and Raginsky (2017), and so are not directly amenable to the single trajectory setting. We also note that information-theoretic generalization bounds have previously found other applications, such as in the analysis of stochastic gradient descent (Neu et al., 2021).

3 Preliminaries and Notation

Results

Our main result is an error bound that controls the distance between the estimate $\hat{f}$, defined by the nonparametric LSE (4), and the ground truth $f_{\star}$ in terms of a fresh sample $\xi$ drawn independently of the algorithm from the mixture distribution over the samples $(x_{0},\dots,x_{T-1})$.

To establish inequality (7), we rely on an information-theoretic decoupling argument given in Proposition 3. Informally, Proposition 3 allows us to decompose the error as

The training error term above is controlled by the offset basic inequality (Lemma 1). Discretizing and proceeding through chaining yields a maximal inequality with $\alpha$ and $\gamma$ as trade-off parameters. The generalization error term above is defined in terms of the discretization parameter $\delta$, which controls the discrepancy between the quantized model $\hat{f}_{\delta}$ (defined in (6)) and the LSE $\hat{f}$. The quantized model $\hat{f}_{\delta}$, used solely in the proof, “generalizes well” since it belongs to a finite hypothesis class by construction.

A similar statement holds for $q\geq 2$ and can be found following the proof of Theorem 2 in Appendix A.1. We also have a version of the above theorem, proven in Appendix A.2, applicable to the parametric entropy growth regime.

This bound shows that more stable systems, as captured by the Lipschitz constant $L_{\star}$, have smaller generalization error. This interpretation is in line with the recent trend of using stability bounds to study learning algorithms applied to data generated by a dynamical system, see for example Boffi et al. (2021) and Tu et al. (2021). Finally we note that although we restricted our analysis to contracting systems, our results are easily extended to a more general notion of nonlinear stability. In particular, we extend Proposition 1 in Appendix B.2 to systems satisfying a notion of exponential incremental input-to-state stability, a standard notion from robust nonlinear control theory (Angeli, 2002).

Mixing systems

2 Summary

Our results show that a large class of nonlinear systems can be learned at the minimax optimal nonparametric rate $T^{-1/(2+q)}$ using the LSE (4). This significantly extends our current understanding of nonasymptotic learning of dynamical systems from single trajectory data. By contrast, previous work assumes i.i.d. data or focuses on either linear models (Simchowitz et al., 2018; Tsiamis and Pappas, 2019; Jedra and Proutiere, 2020) or parametric models with known nonlinearities (Foster et al., 2020; Sattar and Oymak, 2020; Mania et al., 2020; Jain et al., 2021).

Proof Strategy for Theorem 1

Our analysis of the least squares estimator (4) begins with the following information-theoretic decoupling estimate inspired by Russo and Zou (2019) and Xu and Raginsky (2017).

where $Z=(x_{0},\dots,x_{T-1},y_{0},\dots,y_{T-1})$ and where $\xi$ has uniform mixture distribution over the covariates $(x_{0},\dots,x_{T-1})$ of system (1) and is independent of all other randomness.

The proof of the estimate (11) relies on the Donsker-Varadhan variational representation of relative entropy and is given in Appendix C.

We now describe our analysis of the in-sample prediction (or training) error, namely the first term appearing on the right hand side of inequality (12). We start with an inequality due to Liang et al. (2015), which is a variant of the basic inequality of least squares and that is crucial to analyzing the in-sample prediction error for correlated data.

This leads to the maximal inequality (15) of Lemma 3 below.

As observed by Liang et al. (2015), if we had not included the offset term $-\|f(x_{t})\|^{2}$, a naive bound would have yielded $\mathbf{E}\sup_{f\in S}M_{T}(f)\lesssim\sqrt{T\ln|S|}$, penalizing us by a factor $\sqrt{T}$ for the scale of $\sum_{t=0}^{T-1}\langle w_{t},f(x_{t})\rangle$.

Applications of Theorem 1

Suppose that we are in the autoregressive setting (2). If

The factor $T^{-1/3}$ is optimal, and appears even in the i.i.d. case, see for example (Wainwright, 2019, Ch. 13).

The next example revisits the generalized linear models using Theorem 3. This model has recently been analyzed using recursive methods (Foster et al., 2020; Sattar and Oymak, 2020; Jain et al., 2021).

for some $C>0$ and where the Frobenius norm of $A$ is given by $\|A\|_{F}=\sqrt{\operatorname{tr}A^{\top}A}$. This setting is a special case of the autoregressive system (2) with $f_{\star}(\cdot)=\phi(A_{\star}\>\cdot\>)$.

A proof of this claim can be found in Appendix D.3. The dependency on $T,d_{x}$ and $L_{\star}$ matches Theorem 2 of Foster et al. (2020), which is optimal in $T$ and $d_{x}$.

Discussion

We have leveraged recently developed information-theoretic tools (Russo and Zou, 2019; Xu and Raginsky, 2017) to analyze the nonparametric LSE (4) for learning dynamical systems. Our analysis yields, to the best of our knowledge, the first rate-optimal bounds for nonparametric estimation of stable or otherwise mixing nonlinear systems from a single trajectory. In addition, our results are able to capture, as a special case, existing parametric rates in the literature (Foster et al., 2020; Sattar and Oymak, 2020).

Finally, an open problem is to determine for which learning problems the system (1) is required to mix in the single trajectory setting. Most previous works on learning in nonlinear dynamical systems rely on similar mixing time or stability arguments. The cost of this is typically a multiplicative factor in the final bound that degrades as stability is lost (Foster et al., 2020; Sattar and Oymak, 2020; Boffi et al., 2021). In contrast, it is well-known that this dependency can be avoided for learning in linear systems (Lai and Wei, 1982; Simchowitz et al., 2018). Recently Jain et al. (2021) showed under a strong invertibility condition that dependency on the mixing time can also be avoided for the generalized linear model (16). This leaves open the question whether learning without mixing is possible in situations beyond the generalized linear model.

Ingvar Ziemann and Henrik Sandberg are supported by the Swedish Research Council (grant 2016-00861). Nikolai Matni is supported in part by NSF awards CPS-2038873 and CAREER award ECCS-2045834, and a Google Research Scholar award.

Appendix A Proof of Theorem 1 and its Corollaries

We now turn to the proof of Theorem 1. First, we begin by applying Proposition 3 to the discretized estimator $\hat{f}$. Let us begin by bounding the generalization error of the quantized estimator $\hat{f}_{\delta}$ as defined by (6). We may write

It remains to bound the in-sample-prediction error. We have

The main technincal chaining step is given in Lemma 4. Namely, by appealing to Lemma 4 and combining with (17) and (18) we find

since $\mathbf{E}\|\hat{f}(\xi)-f_{\star}(\xi)\|_{2}\leq\mathbf{E}\|\hat{f}_{\delta}(\xi)-f_{\star}(\xi)\|_{2}+\delta$ by the triangle inequality. The result follows after pulling the factor $2\delta^{2}$ out of the square root sign using the triangle inequality. $\blacksquare$

For large spaces and fine grained coverings, the metric entropy starts to dominate the scale free process $M_{T}(f)$ appearing in Lemma 3. The analysis of $M_{T}(f)$ in Lemma 4 below essentially follows that in Liang et al. (2015) (compare with their Lemma 6) with certain slight simplifications due to the added structure the uniform topology on $C(\mathsf{X}\to\mathsf{Y})$ affords us. We begin with an analogue of Lemma 3 which takes the scale of the functions considered into account.

As noted in Liang et al. (2015), the optimal value for $\gamma$ in Lemma 4 is of the same nature as when obtained by other methods, see for example Chapter 13 of Wainwright (2019) for a more standard approach.

Observe first that for any fixed $\alpha>0$ we have, simply by discarding the negative second order term, Cauchy-Schwarz, and a standard subgaussian concentration inequality for $\mathbf{E}\|w_{t}\|_{2}$:

A standard one-step discretization bound (c.f. the proof of Proposition 5.17 in Wainwright (2019)) combined with the finite class maximal inequality of Lemma 3 yields for fixed $\gamma>0$:

Having extracted the fast rate term for scales larger than $\gamma$, we proceed with a chaining bound on the second term above. Since $M_{T}(f)$ satisfies the maximal inequality (23) with $r=\gamma$, chaining (as in Theorem 5.22 of Wainwright (2019)) yields

Under the hypothesis (8) we may use Theorem 1 with $\alpha=0$ to write

The claim follows by solving for the optimal balance

A.2 Proof of Theorem 3

By virtue of Theorem 1 and by selecting $\gamma=\alpha$ we have

Let now $\alpha=1/\sqrt{d_{y}}\sigma_{w}T^{2}$ and $\delta=1/T$. Then (22) by using the hypothesis (9) becomes

A.3 Proof of Auxilliary Results

By optimality of $\hat{f}$ to the prediction error objective we have that

Rearranging and expanding the square gives the basic inequality

which after multiplying both sides by $2$ can be rearranged again to give

Proof of Lemma 3

The proof is a straight-forward modification of the standard proof for bounding the expected supremum of subgaussian maxima. By Jensen’s inequality and monotonicity of the exponential it follows that

Choosing $\lambda=1/2\sigma_{w}^{2}$, application of Lemma 2 yields $\exp\left(\frac{1}{2\sigma^{2}}\mathbf{E}\sup_{f\in S}M_{T}(f)\right)\leq|S|$ which is equivalent to the result. $\blacksquare$

Proof of Lemma 2

Fix $\lambda>0$. By Jensen’s inequality and monotonicity of the exponential it follows that

Using $\|f\|_{\infty}\leq r$ and the tower property let us now estimate

Hence after applying logarithms to both sides of equation (24), we find

which yields the result after optimizing over $\lambda>0$. ∎

A.4 Finite Classes

Appendix B Proofs Related to Stability and Learning

with the convention $\mathbf{E}_{t-1}=\mathbf{E}$. Note now that $F(x_{0},\dots,x_{T-1})-\mathbf{E}F(x_{0},\dots,x_{T-1})=\sum_{t=0}^{T-1}\Delta_{t}$ so to arrive at the desired conclusion we need to prove that the $\Delta_{s}$ are uniformly bounded.

To this end, for a fixed $s$, we define two couplings of $(x_{s},\dots,x_{T-1})$ via

which vary only in their initial condition $z,z^{\prime}$ but are constructed with the same sequence $(w_{s},\dots,w_{T-1})$. We now compute

where the first inequality uses the Markov property to realize the conditional expectations as functions of $x_{s}$ and $x_{s-1}$ respectively. The other inequalities follow by application of the triangle inequality and the $2L$-Lipschitzness of $h$.

Let us now bound the $\|\cdot\|_{\mathsf{X}}$-distance between $z_{t}$ and $z_{t}^{\prime}$:

Combining equations (26) and (27), and noting that a symmetric argument applies to $-\Delta_{s}$ it follows that

Expressing $F$ as a telescoping sum over $\Delta_{s}$, we can compute its moment generating function in combination with the tower property:

using Hoeffding’s inequality to bound the conditional moment generating functions of the bounded random variables $\Delta_{j}$ using the inequality (28) (see Hoeffding (1963) or Example 2.4. in Wainwright (2019)). $\blacksquare$

B.2 Extension to Exponential Incremental Input-to-State Stability

with $\phi_{t},\psi_{t}\in\mathsf{X},\eta_{t},\zeta_{t}\in\mathsf{W}$ it holds for all $t\in[T]$ that

by repeated application of the triangle inequality and since $l$ is $L$-Lipschitz.

Let us now bound the $\rho_{\mathsf{X}}$-distance between $x_{t}$ and $z_{t}$ under the hypothesis that $\eta_{t}=\zeta_{t},\forall t\neq j$. Then we have using the E-$\delta$-ISS bound in equation (29) that

We may proceed with the analysis by defining the martingale difference sequence

which has bounded absolute value by independence of the sequence $\{\eta_{t}\}$ and (33). Observe that this allows us to express $F$ as a telescoping sum, which we can readily use to compute the moment generating function in combination with the tower property:

using Hoeffding’s inequality to bound the conditional moment generating functions of the bounded random variables $\Delta_{j}$ using (33) (see Hoeffding (1963) or Example 2.4. in Wainwright (2019)). ∎

B.3 Proof of Proposition 2

Appendix C Proof of the Decoupling Estimate, Proposition 3

In what follows, we compare probability integrals under different distributions. More precisely, we wish to relate the joint distribution of the least squares estimator (4) and the samples from the system (1) with the product measure of their marginals. The following variational formulation of $D(P\|Q)$, due to Donsker and Varadhan (1975), is key:

Moreover, if $D(\mathbf{P}\|\mathbf{Q})<\infty$, then equality in (34) is attained at $F=\log\frac{d\mathbf{P}}{d\mathbf{Q}}$.

Equipped with Lemma 6, and inspired by the work of Russo and Zou (2019) and Xu and Raginsky (2017), we now turn to the proof of Proposition 3. We remark that the first paragraph of the proof is identical to the proof of Lemma 1 in Xu and Raginsky (2017). As it is central to our argument, we reproduce it below.

We begin by observing that by rescaling $F$ in (34) by $\lambda$, we obtain

For any $F$ which is $\sigma^{2}$-subgaussian under $\mathbf{Q}$, we have that

Combining inequalities (35) and (36), we see that

which after choice of $\lambda=-\frac{\sqrt{D(P\|Q)}}{\sqrt{2\sigma^{2}}}$ and rearranging becomes

We now specialize this known result to our setting. Let us now choose

Observe that for $\mathbf{P},\mathbf{Q}$ as above, $D(\mathbf{P}\|\mathbf{Q})=I((f,g);Z)$. Let further $(x^{\prime}_{0},\dots,x^{\prime}_{T-1})$ be equal in distribution to $(x_{0},\dots,x_{T-1})$ but independent from $f$ and $g$. In other words $(x_{0},\dots,x_{T-1},f,g)$ is drawn from $\mathbf{P}$ and $(x^{\prime}_{0},\dots,x^{\prime}_{T-1},f,g)$ is drawn from $\mathbf{Q}$. Let also $\tau$ be uniformly distributed over $\{0,\dots,T-1\}$ and independent of all other randomness so that we may take $\xi=x^{\prime}_{\tau}$. Hence, for these choices, inequality (37) combined with Jensen’s inequality yield

by linearity of expectation and reformulating the mixture component. Inequality (i) follows from inequality (37) and inequality (ii) from Jensen’s inequality. $\blacksquare$

C.1 Extension: Generalization Bounds for Dynamical Systems

The problem of statistical learning is to find a hypothesis $h\in\mathsf{H}$ that minimizes

with $Z=(x_{0},\dots,x_{T-1},y_{0},\dots,y_{T-1})$ and where $\mathbf{E}_{Z}$ denotes integration over the randomness in $Z$. Let $H$ be a randomized learning algorithm (a random, data-dependent element of $\mathsf{H}$). We define its generalization error by

where $\bar{Z}$ is equal to $Z$ in distribution but independent of $H$. By combining Lemma 1 of Xu and Raginsky (2017) with Proposition 4 we arrive at the following inequality.

Suppose that $l$ is $L$-lipschitz in its first two arguments:

that $(x_{0},\dots,x_{T-1})$ is $(a,b,r)$ E-$\delta$ISS and that $(y_{0},\dots,y_{T-1})$ is $(a^{\prime},b^{\prime},r^{\prime})$ E-$\delta$ISS. Then

where $r_{\max}=\max(r,r^{\prime})$ and $B=\sup_{w,w^{\prime}\in\mathsf{W}}\rho_{\mathsf{W}}(w,w^{\prime})$.

In principle a direct proof using the methods from Appendix B is possible. For brevity, we instead show how the result can be reduced to the statement of Proposition 4.

Let $\mathsf{X}^{\prime}=\mathsf{X}\times\mathsf{Y}$ and define the extended dynamics $\phi^{1}_{t}=x_{t},\phi^{2}_{t}=y_{t}$. Then

or $\phi_{t+1}=G_{t}(\phi_{t},w_{t})$ in brief. Since $G^{x}$ and $G^{y}$ are both E-$\delta$ISS as system from $(\mathsf{W},\rho_{\mathsf{W}})$ to $(\mathsf{X},\rho_{\mathsf{X}})$ and $(\mathsf{Y},\rho_{\mathsf{Y}})$ respectively, it follows that $G$ is $(\max(a,a^{\prime}),2b,\max(r,r^{\prime}))$ E-$\delta$ISS from $(\mathsf{W},\rho_{\mathsf{W}})$ to $(\mathsf{X^{\prime}},\rho_{\mathsf{X^{\prime}}})$ with $\rho_{\mathsf{X}^{\prime}}=\rho_{\mathsf{X}}+\rho_{\mathsf{Y}}$. Hence, we may apply Proposition 4 to conclude that

is $\frac{32L^{2}B^{2}b^{2}}{(1-r_{\max})^{2}T}$-subgaussian for each fixed $h$ where $r_{\max}=\max(r,r^{\prime})$. The result follows by applying Lemma 1 of Xu and Raginsky (2017). ∎

Appendix D Supporting Material for the Examples in Section 4

We may assume that the eigenvalues $\{\lambda_{i}\}$ are ordered as $\lambda_{1}\geq\lambda_{2}\geq\dots$. Fix an integer $m$ and define

D.2 An Experiment Supporting Example 2

We then use $f_{\star}$ to generate training trajectories of varying length to be used in the LSE (4), as well as use $f_{\star}$ to generate $500$ i.i.d. draws from the stationary distributionApproximated by running the system for a burn in time of $T=1000$ time-steps before sampling from it.. To implement the LSE (4) we pass by the dual problem, kernel ridge regression, to estimate $\hat{f}$. We then approximate the $2$-norm distance $\|\hat{f}(\xi)-f_{\star}(\xi)\|_{2}$ by drawing $1000$ fresh trajectories of length and averaging over the final sample. We average our results over $10$ independent systems (random draws of $f_{\star}$) and plot our experiment in Figure 1. It is interesting to note that the slope of the logarithmic plot is slightly less steep than $-1/2$. This is consistent with the near parametric rate of convergence suggested by Example 2 and the exponential eigenvalue decay of the kernel $K(x,z)=\exp\left(-\frac{1}{2}\|x-z\|_{2}^{2}\right)$, see Wainwright (2019), page 399.

D.3 Proof of the claim in Example 3

References