# Learning Recurrent Neural Net Models of Nonlinear Systems

Joshua Hanson University of Illinois; e-mail<EMAIL_ADDRESS>Maxim
Raginsky University of Illinois; e-mail<EMAIL_ADDRESS>Eduardo Sontag
Northeastern University; e-mail<EMAIL_ADDRESS>
###### Abstract

We consider the following learning problem: Given sample pairs of input and
output signals generated by an unknown nonlinear system (which is not assumed
to be causal or time-invariant), we wish to find a continuous-time recurrent
neural net with hyperbolic tangent activation function that approximately
reproduces the underlying i/o behavior with high confidence. Leveraging
earlier work concerned with matching output derivatives up to a given finite
order (Sontag, 1998), we reformulate the learning problem in familiar system-
theoretic language and derive quantitative guarantees on the sup-norm risk of
the learned model in terms of the number of neurons, the sample size, the
number of derivatives being matched, and the regularity properties of the
inputs, the outputs, and the unknown i/o map.

##  1 Introduction

We consider a learning-theoretic framework for system identification, where
the goal is to approximate an unknown nonlinear input/output (i/o) map by a
continuous-time recurrent neural net (RNN) on the basis of a finite collection
of input/output pairs. The approximation criterion is the
$\mathcal{L}^{\infty}$ norm of the difference between the ground-truth output
and the one predicted by the learned model.

Earlier work (Sontag, 1998) has addressed a related problem of reproducing
output $k$-jets as a function of input $(k-1)$-jets via RNNs, where a $k$-jet
is defined as the vector of derivatives up to order $k$ of the output
(respectively, input) signal evaluated at the initialization time.
Equivalently, $k$-jets can be defined as optimal degree-$k$ polynomial
approximations, i.e., truncated Taylor series. In this work, we build on the
aforementioned result about matching $k$-jets, but work in a familiar system-
theoretic setting using the language of i/o maps.

The training procedure used both in that work and here differs to some extent
from the standard “backpropagation through time” algorithm often used for
training RNNs. Rather than forward-propagating training inputs through the net
or setting up an adjoint equation for the net parameters, we “pull back” the
input signals to the initialization time and compute the corresponding output
$k$-jets, which can be expressed explicitly as a function of the weights of
the neural net, its initial state, and the $(k-1)$-jets of the input. The loss
function is then evaluated with the predicted and ground-truth output jets. At
training time, there is no requirement to advance inputs through the network.
In essence, the RNN can be trained using conventional nonlinear regression as
if it were a single-layer feedforward net.

In the following section, we describe the learning problem, state our
assumptions, and introduce the needed notation. Proceeding in Section 3, we
develop the setting of jets and use this language to formulate the proposed
learning procedure and state our main results. Some directions for future
research are also outlined. Section 4.1 contains some technical lemmas,
followed by the proof of the main theorems in Sections 4.2 and 4.3.

##  2 Problem formulation

In this work, a system (or an i/o map) is a nonlinear operator
$\mathsf{F}\mathrel{\mathop{\mathchar 58\relax}}C([0,T])\to C([0,T])$, where
$C([0,T])$ is the Banach space of continuous functions
$u\mathrel{\mathop{\mathchar 58\relax}}[0,T]\to{\mathbb{R}}$ equipped with the
sup norm

$\displaystyle\|u\|_{\infty}\mathrel{\mathop{\mathchar
58\relax}}=\sup_{t\in[0,T]}|u(t)|.$

The learning problem can be phrased as follows: Let $N$ pairs
$(u^{1},y^{1}),\ldots,(u^{N},y^{N})$ be given, where $u^{i}\in C([0,T])$ are
the inputs and $y^{i}=\mathsf{F}u^{i}\in C([0,T])$ are the corresponding
outputs. We wish to construct a system
$\hat{\mathsf{F}}\mathrel{\mathop{\mathchar 58\relax}}C([0,T])\to C([0,T])$
that approximately reproduces the unknown system $\mathsf{F}$ on a given class
of inputs $\mathcal{U}\subset C([0,T])$. (We focus on single-input, single-
output systems mainly for simplicity; our approach easily extends to multiple
inputs and outputs.) Next, we impose a set of minimal assumptions on the
system $\mathsf{F}$ and specify the accuracy criterion.

We first recall the definition of the modulus of continuity of a function
$u\in C([0,T])$ (see, e.g., Ch. 2 of DeVore and Lorentz (1993)):

$\displaystyle\omega_{u}(\delta)\mathrel{\mathop{\mathchar
58\relax}}=\sup_{t_{1},t_{2}\in[0,T]\atop|t_{1}-t_{2}|\leq\delta}|u(t_{1})-u(t_{2})|.$

The function $\delta\mapsto\omega_{u}(\delta)$ is nondecreasing, and, since
any $u\in C([0,T])$ is uniformly continuous, $\lim_{\delta\downarrow
0}\omega_{u}(\delta)=\omega_{u}(0)=0$. (In fact, we will refer to any function
with these properties that majorizes $\omega_{u}$ as a modulus of continuity
of $u$.) We impose the following assumption on the class of inputs
$\mathcal{U}$:

######  Assumption 2.1.

The inputs in $\mathcal{U}$ are uniformly bounded, i.e.,
$R\mathrel{\mathop{\mathchar
58\relax}}=\sup_{u\in\mathcal{U}}\|u\|_{\infty}<\infty$, and equicontinuous
with common modulus of continuity $\omega_{\mathcal{U}}(\delta)$:

$\sup_{u\in\mathcal{U}}\omega_{u}(\delta)\leq\omega_{\mathcal{U}}(\delta),$

where $\omega_{\mathcal{U}}(\delta)$ is a nondecreasing function that
satisfies $\lim_{\delta\downarrow
0}\omega_{\mathcal{U}}(\delta)=\omega_{\mathcal{U}}(0)=0$.

We also use moduli of continuity to describe the regularity of $\mathsf{F}$:

######  Assumption 2.2.

The output space $\mathcal{Y}=\mathsf{F}(\mathcal{U})$, i.e., the image of
$\mathcal{U}$ under $\mathsf{F}$, is equicontinuous with a common modulus of
continuity $\omega_{\mathcal{Y}}(\delta)$.

Finally, we assume that $\mathsf{F}$ has the bounded-in, bounded-out property:

######  Assumption 2.3.

$\gamma_{\mathsf{F}}(R)\mathrel{\mathop{\mathchar
58\relax}}=\displaystyle\sup_{u\in C[0,T]\atop\|u\|_{\infty}\leq
R}\|\mathsf{F}u\|_{\infty}<\infty$.

Examples of classes of inputs that satisfy Assumption 2.1 include:

  * •

finite combinations of Fourier terms

$\displaystyle u(t)=\sum^{m}_{i=1}c_{i}\sin(\omega_{i}t+\alpha_{i}),$ (2.1)

for all $m\geq 1$, provided the coefficients $c_{i}$ and the frequencies
$\omega_{i}$ satisfy the inequalities $\sum^{m}_{i=1}|c_{i}|\leq R$ and
$\sum^{m}_{i=1}|c_{i}\omega_{i}|\leq L$ for some fixed finite constants $R$
and $L$;

  * •

polynomial inputs

$\displaystyle u(t)=\sum^{m}_{i=0}c_{i}t^{i}$ (2.2)

for all $m\geq 0$, provided the coefficients $c_{i}$ satisfy the inequalities
$\sum^{m}_{i=0}|c_{i}|T^{i}\leq R$ and $\sum^{m}_{i=1}i|c_{i}|T^{i-1}\leq L$
for some $R,L<\infty$.

In both cases, Assumption 2.1 holds with
$\omega_{\mathcal{U}}(\delta)=L\delta$. As an example of a system satisfying
our hypotheses, consider a state-space model of the form

$\displaystyle\dot{x}(t)$ $\displaystyle=f(x(t))+g(x(t))u(t),$ $\displaystyle
y(t)$ $\displaystyle=h(x(t))$

where $x(t)\in{\mathbb{R}}^{n}$ is an internal state with a given initial
condition $x(0)=\xi$. Let $\mathsf{F}$ be the induced i/o map that sends the
input $u\mathrel{\mathop{\mathchar 58\relax}}[0,T]\to{\mathbb{R}}$ to the
output $y\mathrel{\mathop{\mathchar 58\relax}}[0,T]\to{\mathbb{R}}$. Then it
is not hard to verify that Assumptions 2.2 and 2.3 will hold if the functions
$f\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$,
$g\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$,
and $h\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}^{n}\to{\mathbb{R}}$
are all Lipschitz continuous.

We assume the sample inputs $u^{i}$ are independent and identically
distributed (i.i.d.) random elements of $C([0,T])$ drawn according to a fixed
Borel probability measure $\mu$ supported on $\mathcal{U}$. Thus, the input-
output pairs $(u^{1},y^{1}),\ldots,(u^{N},y^{N})$, with
$y^{i}=\mathsf{F}u^{i}$, are themselves i.i.d. random elements of the product
space $\mathcal{U}\times\mathcal{Y}$. For instance, $\mathsf{F}$ could be a
model of a two-terminal electronic device whose terminals are connected
through a switch to an excitation circuit consisting of linear and nonlinear
elements and current and/or voltage sources (Chua, 1980). Suppose this
excitation circuit is specified by a vector of parameters
$\theta\in{\mathbb{R}}^{d}$; e.g., it could be used to generate sinusoidal or
polynomial inputs, as in (2.1) or (2.2). We can then generate $N$ i.i.d.
samples $\theta^{1},\ldots,\theta^{N}$ according to a fixed Borel probability
measure $\nu$ on ${\mathbb{R}}^{d}$. Each sample $\theta^{i}$ corresponds to a
random realization of the excitation circuit so that, when the switch is
closed at time $t=0$ and open at time $t=T$, we can take the input
$u^{i}\mathrel{\mathop{\mathchar 58\relax}}[0,T]\to{\mathbb{R}}$ to be the
voltage waveform across $\mathsf{F}$ and the output
$y^{i}\mathrel{\mathop{\mathchar 58\relax}}[0,T]\to{\mathbb{R}}$ to be the
current waveform through $\mathsf{F}$. (This assumes that $\mathsf{F}$ is
voltage-controlled.) Thus, the probability measure on the input-output space
$\mathcal{U}\times\mathcal{Y}$ is well-defined but specified indirectly
through $\nu$, the structure of the excitation circuit, and $\mathsf{F}$. At
any rate, given $\mathsf{F}$ and an approximating system $\hat{\mathsf{F}}$,
we define the $\mathcal{L}^{\infty}$ risk

$\displaystyle\mathcal{L}(\hat{\mathsf{F}})\mathrel{\mathop{\mathchar
58\relax}}={\mathbf{E}}_{\mu}\big{[}\|\hat{\mathsf{F}}u-\mathsf{F}u\|_{\infty}\big{]},$

where the expectation is taken with respect to the probability measure $\mu$
on the input space $\mathcal{U}$. The goal is to generate, on the basis of the
observed input-output pairs $(u^{i},y^{i})$, an approximate system
$\hat{\mathsf{F}}$ from a given model class $\mathcal{F}$, so that the risk
$\mathcal{L}(\hat{\mathsf{F}})$ (which is a random variable due to its
dependence on the training data) is small with high probability.

##  3 The proposed learning procedure and its performance

#### Recurrent neural nets.

We first specify the model class $\mathcal{F}$ from which our learning
procedure will select the approximation $\hat{\mathsf{F}}$. The systems in
$\mathcal{F}$ are described by differential equations of the form

$\displaystyle\dot{x}(t)$ $\displaystyle=\sigma^{(n)}(Ax(t)+bu(t))$ (3.1a)
$\displaystyle y(t)$ $\displaystyle=c^{\hbox{{\tiny T}}}x(t)$ (3.1b)

for $t\in[0,T]$ with initial condition $x(0)=\xi\in{\mathbb{R}}^{n}$. Here,
$A\in{\mathbb{R}}^{n\times n}$, $b,c\in{\mathbb{R}}^{n}$, and
$\sigma^{(n)}\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ is the diagonal map
$\sigma^{(n)}((x_{1},\ldots,x_{n})^{\hbox{{\tiny
T}}})\mathrel{\mathop{\mathchar
58\relax}}=(\sigma(x_{1}),\ldots,\sigma(x_{n}))^{\hbox{{\tiny T}}}$ with
$\sigma(r)=\tanh(r)$. The system (3.1) is a continuous-time recurrent neural
net (RNN) with $n$ neurons. Each pair $(\Sigma,\xi)$, where
$\Sigma\mathrel{\mathop{\mathchar 58\relax}}=(A,b,c)$, specifies an input-
output map $\mathsf{F}_{\Sigma,\xi}$ that sends an input $u\in C([0,T])$ to
the corresponding output $y$ given by

$\displaystyle y(t)=c^{\hbox{{\tiny T}}}\xi+\int^{t}_{0}c^{\hbox{{\tiny
T}}}\sigma^{(n)}(Ax(\tau)+Bu(\tau))\operatorname{d\\!}\tau,\qquad t\in[0,T].$

For $0<M<\infty$, we define the class of systems

$\displaystyle\mathcal{F}(M)\mathrel{\mathop{\mathchar
58\relax}}=\left\\{\mathsf{F}_{\Sigma,\xi}=\mathsf{F}_{(A,b,c),\xi}\mathrel{\mathop{\mathchar
58\relax}}\|A\|,|b|,|c|,|\xi|\leq M\right\\},$

where $\|A\|$ is the spectral norm of $A$, and $|b|,|c|,|\xi|$ are the
$\ell^{2}$ norms of $b,c,\xi$. Our learning procedure will generate a pair
$(\hat{\Sigma},\hat{\xi})$ based on the data $\\{(u^{i},y^{i})\\}$, and output
the model
$\hat{\mathsf{F}}=\mathsf{F}_{\hat{\Sigma},\hat{\xi}}\in\mathcal{F}(M)$. With
a slight abuse of notation, we will often write
$(\Sigma,\xi)\in\mathcal{F}(M)$ instead of
$\mathsf{F}_{\Sigma,\xi}\in\mathcal{F}(M)$.

#### Jets.

While it is well-known that RNNs of the form (3.1) are universal approximators
for i/o maps $\mathsf{F}$ that admit smooth nonlinear state-space realizations
$\dot{x}=f(x,u)$, $y=g(x)$ (Sontag, 1992; Funahashi and Nakamura, 1993; Hanson
and Raginsky, 2020), here we are not assuming that $\mathsf{F}$ admits such a
realization (in fact, we are not even requiring $\mathsf{F}$ to be causal or
time-invariant). Nevertheless, we will show that, provided Assumptions 2.1–2.3
hold for our learning problem, we will be able to approximate $\mathsf{F}$ by
a recurrent net model which will have the properties of causality and time
invariance by construction. Our approach proceeds by way of reducing the
infinite-dimensional problem of learning the i/o map
$\mathsf{F}_{\hat{\Sigma},\hat{\xi}}$ to a certain finite-dimensional problem
(Sontag, 1998). To that end, consider a system of the form (3.1) fed with an
input $u$ which has at least $k-1$ derivatives at $t=0$. Then the output
$y=\mathsf{F}_{\Sigma,\xi}u$ will have at least $k$ derivatives at $t=0$,
which can be computed explicitly as

$\displaystyle y^{(0)}(0)$ $\displaystyle=c^{\hbox{{\tiny T}}}\xi,\quad
y^{(1)}(0)=c^{\hbox{{\tiny T}}}\sigma^{(n)}(A\xi+Bu(0)),\quad\ldots.$

We can then define the map $Y_{k,\Sigma,\xi}\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{k}\to{\mathbb{R}}^{k+1}$ according to

$\displaystyle
Y_{k,\Sigma,\xi}((u(0),u^{\prime}(0),\ldots,u^{(k-1)}(0))^{\hbox{{\tiny
T}}})\mathrel{\mathop{\mathchar
58\relax}}=(y(0),y^{\prime}(0),\ldots,y^{(k)}(0))^{\hbox{{\tiny T}}},$

where
$y^{(i)}(0)=\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\big{|}_{t=0}\mathsf{F}_{\Sigma,\xi}u(t)$
for $0\leq\ell\leq k$. We can also phrase this in terms of jets, where the
$k$-jet at $t=0$ of a function $f\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}\to{\mathbb{R}}$ which is $C^{k}$ in some neighborhood
of $t=0$ is the degree-$k$ polynomial

$\displaystyle J^{k}_{0}f(s)\mathrel{\mathop{\mathchar
58\relax}}=\sum^{k}_{\ell=0}f^{(\ell)}(0)\frac{s^{\ell}}{\ell!}.$

Then, for inputs that are $C^{k-1}$ in some neighborhood of $t=0$, the map
$Y_{k,\Sigma,\xi}\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{k}\to{\mathbb{R}}^{k+1}$ can be lifted to a map from
${(k-1)}$-jets to $k$-jets via $J^{k-1}_{0}u\mapsto
J^{k}_{0}\mathsf{F}_{\Sigma,\xi}u$, where the vector of coefficients of
$J^{k}_{0}\mathsf{F}_{\Sigma,\xi}u$ is the image of the vector of coefficients
of $J^{k-1}_{0}u$ under $Y_{k,\Sigma,\xi}$. Thus, at least for inputs $u$ that
are of class $C^{k-1}$, the i/o map $\mathsf{F}_{\Sigma,\xi}$ will be a good
approximation to $\mathsf{F}$ provided (a) the outputs $\mathsf{F}u$ and
$\mathsf{F}_{\Sigma,\xi}u$ can be accurately approximated by their $k$-jets
$J^{k}_{0}\mathsf{F}u$ and $J^{k}_{0}\mathsf{F}_{\Sigma,\xi}u$ (e.g., if $k$
is sufficiently large) and (b) the $k$-jets $J^{k}_{0}\mathsf{F}u$ and
$J^{k}_{0}\mathsf{F}_{\Sigma,\xi}u$ are close to one another in sup norm on
$[0,T]$.

#### The learning procedure.

Since the inputs in $\mathcal{U}$ and the corresponding outputs in
$\mathcal{Y}$ are only assumed to be continuous, we will first approximate
them by functions that have as many derivatives at $t=0$ as needed. To that
end, we start by defining for each $k\in{\mathbb{N}}$ two linear maps
$\mathsf{S}_{k}\mathrel{\mathop{\mathchar
58\relax}}C([0,T])\to{\mathbb{R}}^{k}$ and
$\mathsf{S}^{*}_{k}\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}^{k}\to
C([0,T])$ according to

$\displaystyle\mathsf{S}_{k}(u)\mathrel{\mathop{\mathchar
58\relax}}=\Big{(}B_{k-1}(u,0),\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}{\Big{|}}_{t=0}B_{k-1}(u,t),\dots,\frac{\operatorname{d\\!}^{\,k-1}}{\operatorname{d\\!}t^{k-1}}{\Big{|}}_{t=0}B_{k-1}(u,t)\Big{)}^{\hbox{{\tiny
T}}},$

and

$\displaystyle\mathsf{S}^{*}_{k}((a_{0},\dots,a_{k-1})^{\hbox{{\tiny
T}}})(t)\mathrel{\mathop{\mathchar
58\relax}}=\sum_{\ell=0}^{k-1}a_{\ell}\frac{t^{\ell}}{\ell!},$

where

$\displaystyle B_{m}(u,t)\mathrel{\mathop{\mathchar
58\relax}}=\sum^{m}_{i=0}u\left(\frac{iT}{m}\right){m\choose
i}\left(\frac{t}{T}\right)^{i}\left(1-\frac{t}{T}\right)^{m-i},\qquad
t\in[0,T]$

is the degree-$m$ Bernstein polynomial of $u\in C([0,T])$ (DeVore and Lorentz,
1993, Ch. 1). Whenever no confusion will arise, we will also write $B_{m}u$
instead of $B_{m}(u,\cdot)$. Observe that
$\mathsf{S}_{k}\circ\mathsf{S}^{*}_{k}=\operatorname{id}_{{\mathbb{R}}^{k}}$
and $(\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u=B_{k-1}(u,\cdot)$. To motivate
the introduction of these objects, we give the following bound on the expected
risk of any recurrent net model $\mathsf{F}_{\Sigma,\xi}$:

######  Theorem 3.1.

Under Assumptions 2.1 and 2.2, the following holds for any $\Sigma=(A,b,c)$,
$\xi$, and $k\geq 2$:

$\displaystyle\begin{split}\mathcal{L}(\mathsf{F}_{\Sigma,\xi})&\leq
2\omega_{\mathcal{Y}}\left(\frac{T}{\sqrt{k}}\right)+2|c||b|e^{\|A\|T}\omega_{\mathcal{U}}\left(\frac{2T}{\sqrt{k}}\right)+|c|Te^{\|A\|T}\sqrt{\frac{n}{k}}\\\
&\qquad\qquad+{\mathbf{E}}_{\mu}\left[\|(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi}\circ\mathsf{S}_{k})u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F})u\|_{\infty}\right].\end{split}$
(3.2)

Note that the last term on the right-hand side of (3.2) is the expectation,
with respect to $u\sim\mu$, of the sup norm of the difference between the
degree-$k$ Bernstein polynomial $B_{k}\mathsf{F}u$ and the degree-$k$
Bernstein polynomial $B_{k}\mathsf{F}_{\Sigma,\xi}B_{k-1}u$.

We are now ready to present our learning procedure. For $1\leq j\leq k$, let
$t_{j}\mathrel{\mathop{\mathchar 58\relax}}=jT/k$, and consider the following
Empirical Risk Minimization (ERM) scheme:

$\displaystyle(\hat{\Sigma},\hat{\xi})\in\operatorname*{arg\,min}_{(\Sigma,\xi)\in\mathcal{F}(M)}\frac{1}{N}\sum^{N}_{i=1}\max_{1\leq
j\leq k}\left|(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi}\circ\mathsf{S}_{k})u(t_{j})-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1})y(t_{j})\right|.$
(3.3)

The objective being minimized in (3.3) is simply the empirical expectation of
the maximum absolute difference between the Bernstein polynomials
$B_{k}\mathsf{F}u$ and $B_{k}\mathsf{F}_{\Sigma,\xi}B_{k-1}u$ on the finite
grid $\\{T/k,2T/k,\ldots,(k-1)T/k,T\\}\subset[0,T]$.

######  Theorem 3.2.

Suppose Assumptions 2.1–2.3 are satisfied, and $N\geq
k(6n^{6}+10n^{3}\log_{2}k)$. Then with probability at least $1-\delta$,

$\displaystyle\begin{split}\mathcal{L}(\hat{\mathsf{F}})&\leq
4\omega_{\mathcal{Y}}\left(\frac{T}{\sqrt{k}}\right)+2M^{2}e^{MT}\omega_{\mathcal{U}}\left(\frac{2T}{\sqrt{k}}\right)+3MTe^{MT}\sqrt{\frac{n}{k}}\\\
&\qquad\qquad+\bar{\mathcal{L}}^{*}+c\left(M(M+\sqrt{n}T)+\gamma_{\mathsf{F}}(R)\right)\sqrt{\frac{k(n^{6}+n^{3}\log_{2}k)\log
N+\log(\frac{1}{\delta})}{N}},\end{split}$ (3.4)

where $c>0$ is an absolute constant and

$\displaystyle\bar{\mathcal{L}}^{*}\mathrel{\mathop{\mathchar
58\relax}}=\inf_{(\Sigma,\xi)\in\mathcal{F}(M)}{\mathbf{E}}_{\mu}\left[\max_{1\leq
j\leq k}\left|(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi}\circ\mathsf{S}_{k})u(t_{j})-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1})y(t_{j})\right|\right].$
(3.5)

The first three terms in (3.4) can be made arbitrarily small by choosing $k$
sufficiently large. The exponential dependence on $M$ and $T$ is an
unavoidable artifact of the Grönwall lemma, and can be removed under
appropriate stability assumptions (Hanson and Raginsky, 2020). The remaining
two terms are, respectively, the approximation error and the estimation error
of the ERM procedure (3.3). While the latter can be made arbitrarily small by
increasing the size $N$ of the training set, the former is an intrinsic
measure of the ability of recurrent neural nets to approximate output $k$-jets
for a randomly chosen input. This minimal error value can be decreased by
considering larger, more expressive nets; a quantitative analysis of this term
is a promising direction for further work. Note also that, for a fixed $k$, we
only need to collect input and output samples on an equispaced grid
$\\{0,T/k,2T/k,\ldots,T\\}$. On the other hand, in many applications (e.g.,
medical or electronic system modeling), the input/output data are often
available at non-uniformly spaced times $0\leq t_{1}<t_{2}<\ldots<t_{k}\leq T$
that are not under the learner’s control. Extending the approach of this paper
to non-uniform (or even random) sampling of time instants is another
interesting future direction.

##  4 Proofs

###  4.1 Technical lemmas

In this section, we collect a few results that will be used in the proofs of
Theorems 3.1 and 3.2.

######  Lemma 4.1.

For any $u\in C([0,T])$ with modulus of continuity $\omega$ and any
$k\in{\mathbb{N}}$, the Bernstein polynomial $B_{k}(u,\cdot)$ has modulus of
continuity $2\omega$ and satisfies

$\displaystyle\|u-B_{k}(u,\cdot)\|_{\infty}\leq
2\omega\left(\frac{T}{\sqrt{k}}\right).$ (4.1)

###### Proof.

The statement about the modulus of continuity of $B_{k}(u,\cdot)$ is a result
of Li (2000). The error bound (4.1) is well-known, but we give a self-
contained probabilistic proof. Let $\bar{\omega}$ be the smallest concave
majorant of $\omega$. Then $\bar{\omega}$ is also a modulus of continuity of
$u$, which further satisfies the inequality
$\omega(\delta)\leq\bar{\omega}(\delta)\leq 2\omega(\delta)$ (DeVore and
Lorentz, 1993, Lemma 6.1).

For any $t\in[0,T]$, we can express the value $B_{k}(u,t)$ as an expectation

$\displaystyle
B_{k}(u,t)={\mathbf{E}}\Big{[}u\Big{(}\frac{T}{k}X_{t}\Big{)}\Big{]}$

where $X_{t}\sim\mathrm{Bin}(k,t/T)$. Then we have the following chain of
(in)equalities:

$\displaystyle|u(t)-B_{k}(u,t)|$
$\displaystyle=\left|u(t)-{\mathbf{E}}\Big{[}u\Big{(}\frac{T}{k}X_{t}\Big{)}\Big{]}\right|$
$\displaystyle\stackrel{{\scriptstyle\mathrm{(a)}}}{{\leq}}{\mathbf{E}}\left[\left|u(t)-u\Big{(}\frac{T}{k}X_{t}\Big{)}\right|\right]$
$\displaystyle\stackrel{{\scriptstyle\mathrm{(b)}}}{{\leq}}{\mathbf{E}}\left[\bar{\omega}\left(\left|\frac{T}{k}X_{t}-t\right|\right)\right]$
$\displaystyle\stackrel{{\scriptstyle\mathrm{(c)}}}{{\leq}}\bar{\omega}\left({\mathbf{E}}\left|\frac{T}{k}X_{t}-t\right|\right)$
$\displaystyle\stackrel{{\scriptstyle\mathrm{(d)}}}{{\leq}}2\omega\left(\frac{T}{\sqrt{k}}\right),$

where (a) is by Jensen’s inequality, (b) follows from the fact that
$\bar{\omega}$ is a modulus of continuity of $u$, (c) is by Jensen’s
inequality and by concavity of $\bar{\omega}$, and (d) uses the monotonicity
of $\omega$ and $\bar{\omega}$, the inequality $\bar{\omega}\leq 2\omega$, and
the inequality ${\mathbf{E}}|U-{\mathbf{E}}U|\leq\sqrt{kp(1-p)}$ for the mean
absolute deviation of a $\mathrm{Bin}(k,p)$ random variable $U$. ∎

######  Lemma 4.2.

Let $\mathsf{G}\mathrel{\mathop{\mathchar 58\relax}}C([0,T])\to C([0,T])$ be a
Lipschitz-continuous i/o map, i.e.,

$\displaystyle\|\mathsf{G}\|_{\mathrm{Lip}}\mathrel{\mathop{\mathchar
58\relax}}=\sup_{u_{1},u_{2}\in C([0,T])\atop u_{1}\neq
u_{2}}\frac{\|\mathsf{G}u_{1}-\mathsf{G}u_{2}\|_{\infty}}{\|u_{1}-u_{2}\|_{\infty}}<\infty.$

Then, for any integer $k\geq 2$ and any $u\in C([0,T])$,

$\displaystyle\|\mathsf{G}u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{G}\circ\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u\|_{\infty}\leq
2\|\mathsf{G}\|_{\mathrm{Lip}}\omega_{u}\left(\frac{2T}{\sqrt{k}}\right)+2\omega_{\mathsf{G}B_{k-1}u}\left(\frac{T}{\sqrt{k}}\right).$

###### Proof.

By the triangle inequality,

$\displaystyle\begin{split}&\phantom{{}={}}\|\mathsf{G}u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{G}\circ\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u\|_{\infty}\\\
&\qquad\leq\|\mathsf{G}u-(\mathsf{G}\circ\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u\|_{\infty}+\|(\mathsf{G}\circ\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{G}\circ\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u\|_{\infty}\\\
&\qquad=\mathrel{\mathop{\mathchar 58\relax}}T_{1}+T_{2}.\end{split}$

Since $(\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u=B_{k-1}u$, Lemma 4.1 gives

$\displaystyle T_{1}$
$\displaystyle\leq\|\mathsf{G}\|_{\mathrm{Lip}}\|u-B_{k-1}u\|_{\infty}\leq
2\|\mathsf{G}\|_{\mathrm{Lip}}\omega_{u}\left(\frac{T}{\sqrt{k-1}}\right)\leq
2\|\mathsf{G}\|_{\mathrm{Lip}}\omega_{u}\left(\frac{2T}{\sqrt{k}}\right),$

where the last inequality follows from the fact that
$\sqrt{k-1}\geq\frac{1}{2}\sqrt{k}$ for $k\geq 2$. Similarly,

$\displaystyle T_{2}$
$\displaystyle=\|\mathsf{G}B_{k-1}u-B_{k}\mathsf{G}B_{k-1}u\|_{\infty}\leq
2\omega_{\mathsf{G}B_{k-1}u}\left(\frac{T}{\sqrt{k}}\right).$

∎

######  Lemma 4.3.

Let $\mathsf{G}$ be the i/o map of the differential dynamical system

$\displaystyle\dot{x}(t)=f(x(t),u(t)),\quad y(t)$
$\displaystyle=h^{\hbox{{\tiny T}}}x(t);\qquad t\in[0,T],\,x(0)=\xi$ (4.2)

with input $u(t)\in{\mathbb{R}}$, state $x(t)\in{\mathbb{R}}^{n}$, and output
$y(t)\in{\mathbb{R}}$. Suppose that $f(\cdot,\cdot)$ is bounded, i.e.,
$\sup_{(x,u)\in{\mathbb{R}}^{n}\times{\mathbb{R}}}|f(x,u)|<\infty$, and
Lipschitz-continuous, i.e.,

$\displaystyle|f(x_{1},u_{1})-f(x_{2},u_{2})|\leq
L_{X}|x_{1}-x_{2}|+L_{U}|u_{1}-u_{2}|.$

Then, for any input $u\in C([0,T])$, the output $y=\mathsf{G}u$ is bounded
with

$\displaystyle\|y\|_{\infty}\leq|h|\left(|\xi|+T\sup_{x,u}|f(x,u)|\right),$
(4.3)

and has modulus of continuity

$\displaystyle\omega_{y}(\delta)\leq\frac{|h|}{L_{X}}\sup_{x,u}|f(x,u)|\left(e^{L_{X}\delta}-1\right).$
(4.4)

Moreover, the i/o map $\mathsf{G}$ is Lipschitz-continuous, with

$\displaystyle\|\mathsf{G}\|_{\mathrm{Lip}}\leq\frac{|h|L_{U}}{L_{X}}\left(e^{L_{X}T}-1\right).$
(4.5)

###### Proof.

For a given input $u\in C([0,T])$ and vector $\xi\in{\mathbb{R}}^{n}$, let
$\varphi^{u}_{s,t}(\xi)$, $0\leq s\leq t\leq T$, denote the flow of (4.2),
i.e., the solution of the ODE

$\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\varphi^{u}_{s,t}(\xi)=f(\varphi^{u}_{s,t}(\xi),u(t)),\qquad
s\leq t\leq T,\,\varphi^{u}_{s,s}(\xi)=\xi.$

The flow map evidently has the semigroup property

$\displaystyle\varphi^{u}_{s,t}(\xi)=\varphi^{u}_{r,t}(\varphi^{u}_{s,r}(\xi)),\qquad
0\leq s\leq r\leq t\leq T.$ (4.6)

In terms of the flow, the i/o map $\mathsf{G}$ is given by
$y(t)=(\mathsf{G}u)(t)=h^{\hbox{{\tiny T}}}\varphi^{u}_{0,t}(\xi)$. To show
that $y$ is bounded, we first estimate

$\displaystyle|\varphi^{u}_{0,t}(\xi)|=\left|\xi+\int^{t}_{0}f(\varphi^{u}_{0,\tau}(\xi),u(\tau))\operatorname{d\\!}\tau\right|\leq|\xi|+T\sup_{x,u}|f(x,u)|,\qquad
0\leq t\leq T$

and then $|y(t)|\leq|h||\varphi^{u}_{0,t}(\xi)|$ by Cauchy–Schwarz. This gives
(4.3).

To obtain the modulus of continuity of $y$, consider two times $0\leq
t_{1}\leq t_{2}\leq T$. Then, using (4.6) and the time invariance of (4.2), we
can write

$\displaystyle|\varphi_{0,t_{2}}^{u}(\xi)-\varphi_{0,t_{1}}^{u}(\xi)|$
$\displaystyle=|\varphi_{t_{1},t_{2}}^{u}(\varphi_{0,t_{1}}^{u}(\xi))-\varphi_{0,t_{1}}^{u}(\xi)|$
$\displaystyle=|\varphi_{0,t_{2}-t_{1}}^{\bar{u}}(\xi^{\prime})-\xi^{\prime}|,$

where $\xi^{\prime}\mathrel{\mathop{\mathchar
58\relax}}=\varphi_{0,t_{1}}^{u}(\xi)$ and $\bar{u}(t)$ is any continuous
extension of the map $[0,t_{2}-t_{1}]\ni t\mapsto u(t+t_{1})$ to $[0,T]$.
Applying the fundamental theorem of calculus and the triangle inequality, we
have

$\displaystyle\phantom{{}={}}|\varphi_{0,t_{2}-t_{1}}^{\bar{u}}(\xi^{\prime})-\xi^{\prime}|$
$\displaystyle\qquad=\Big{|}\int_{0}^{t_{2}-t_{1}}f(\varphi_{0,\tau}^{\bar{u}}(\xi^{\prime}),\bar{u}(\tau))\operatorname{d\\!}\tau\Big{|}$
$\displaystyle\qquad\leq\int_{0}^{t_{2}-t_{1}}\big{(}|f(\varphi_{0,\tau}^{\bar{u}}(\xi^{\prime}),\bar{u}(\tau))-f(\xi^{\prime},\bar{u}(\tau))|+|f(\xi^{\prime},\bar{u}(\tau))|\big{)}\operatorname{d\\!}\tau$
$\displaystyle\qquad\leq\int_{0}^{t_{2}-t_{1}}\big{(}L_{X}|\varphi_{0,\tau}^{\bar{u}}(\xi^{\prime})-\xi^{\prime}|+\sup_{x,u}|f(x,u)|\big{)}\operatorname{d\\!}\tau$
$\displaystyle\qquad=L_{X}\int_{0}^{t_{2}-t_{1}}\big{(}|\varphi_{0,\tau}^{\bar{u}}(\xi^{\prime})-\xi^{\prime}|+\frac{1}{L_{X}}\sup_{x,u}|f(x,u)|\big{)}\operatorname{d\\!}\tau.$

Grönwall’s inequality (Hirsch and Smale, 1974) applied to the function

$\mu(t)\mathrel{\mathop{\mathchar
58\relax}}=|\varphi_{0,t}^{\bar{u}}(\xi^{\prime})-\xi^{\prime}|+\frac{1}{L_{X}}\sup_{x,u}|f(x,u)|$

yields

$\mu(t)\leq\sup_{x,u}|f(x,u)|\frac{1}{L_{X}}e^{L_{X}(t_{2}-t_{1})}.$

Hence, we have

$\displaystyle|y(t_{1})-y(t_{2})|$ $\displaystyle=|h^{\hbox{{\tiny
T}}}(\varphi^{u}_{0,t_{1}}(\xi)-\varphi^{u}_{0,t_{2}}(\xi))|$
$\displaystyle\leq|h||\varphi^{u}_{0,t_{1}}(\xi)-\varphi^{u}_{0,t_{2}}(\xi)|$
$\displaystyle\leq|h|\sup_{x,u}|f(x,u)|\frac{1}{L_{X}}\big{(}e^{L_{X}(t_{2}-t_{1})}-1\big{)}.$

Interchanging the roles of $t_{1}$ and $t_{2}$, we obtain (4.4).

For any time $t\in[0,T]$, the difference in flows due to distinct inputs
$u_{1},u_{2}\in C([0,T])$ is bounded by

$\displaystyle\phantom{{}={}}|\varphi_{0,t}^{u_{1}}(\xi)-\varphi_{0,t}^{u_{2}}(\xi)|$
$\displaystyle\qquad=\Big{|}\int_{0}^{t}f(\varphi_{0,\tau}^{u_{1}}(\xi),u_{1}(\tau))-f(\varphi_{0,\tau}^{u_{2}}(\xi),u_{2}(\tau))\operatorname{d\\!}\tau\Big{|}$
$\displaystyle\qquad\leq\int_{0}^{t}\big{(}|f(\varphi_{0,\tau}^{u_{1}}(\xi),u_{1}(\tau))-f(\varphi_{0,\tau}^{u_{2}}(\xi),u_{1}(\tau))|+|f(\varphi_{0,\tau}^{u_{2}}(\xi),u_{1}(\tau))-f(\varphi_{0,\tau}^{u_{2}}(\xi),u_{2}(\tau))|\big{)}\operatorname{d\\!}\tau$
$\displaystyle\qquad\leq\int_{0}^{t}\big{(}L_{X}|\varphi_{0,t}^{u_{1}}(\xi)-\varphi_{0,t}^{u_{2}}(\xi)|+L_{U}|u_{1}(\tau)-u_{2}(\tau))|\big{)}\operatorname{d\\!}\tau$
$\displaystyle\qquad\leq
L_{X}\int_{0}^{t}\big{(}|\varphi_{0,t}^{u_{1}}(\xi)-\varphi_{0,t}^{u_{2}}(\xi)|+\frac{L_{U}}{L_{X}}\|u_{1}-u_{2}\|_{\infty}\big{)}\operatorname{d\\!}\tau.$

Grönwall’s inequality applied to the function

$\mu(t)\mathrel{\mathop{\mathchar
58\relax}}=|\varphi_{0,t}^{u_{1}}(\xi)-\varphi_{0,t}^{u_{2}}(\xi)|+\frac{L_{U}}{L_{X}}\|u_{1}-u_{2}\|_{\infty}$

yields

$\mu(t)\leq\frac{L_{U}}{L_{X}}e^{L_{X}t}\|u_{1}-u_{2}\|_{\infty},$

hence we have

$|\varphi_{0,t}^{u_{1}}(\xi)-\varphi_{0,t}^{u_{2}}(\xi)|\leq\frac{L_{U}}{L_{X}}\big{(}e^{L_{X}t}-1\big{)}\|u_{1}-u_{2}\|_{\infty}.$

Thus,

$\displaystyle\|\mathsf{G}u_{1}-\mathsf{G}u_{2}\|_{\infty}$
$\displaystyle=\sup_{t\in[0,T]}\left|h^{\hbox{{\tiny
T}}}(\varphi_{0,t}^{u_{1}}(\xi)-\varphi_{0,t}^{u_{2}}(\xi))\right|$
$\displaystyle\leq\frac{|h|L_{U}}{L_{X}}\big{(}e^{L_{X}T}-1\big{)}\|u_{1}-u_{2}\|_{\infty},$

and (4.5) follows. ∎

We also need the following result on the Rademacher averages of VC-subgraph
classes (Farrell et al., 2020). Recall that a class $\mathcal{G}$ of
measurable functions $g\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{d}\to{\mathbb{R}}$ is a VC-subgraph class (Giné and
Nickl, 2016, Sec. 3.6.2) if the class of all sets of the form
$\\{(x,r)\in{\mathbb{R}}^{d}\times{\mathbb{R}}\mathrel{\mathop{\mathchar
58\relax}}g(x)\geq r\\}$ with $g\in\mathcal{G}$ is a Vapnik–Chervonenkis (or
VC) class, i.e., there exists a finite $D\in{\mathbb{N}}$, such that, for each
$m\leq D$, there exist $m$ points $(x^{1},r^{1}),\ldots,(x^{m},r^{m})$ that
are shattered by $\mathcal{G}$, i.e.,

$\displaystyle\\{{\mathbf{1}}_{\\{g(x^{1})\geq
r^{1}\\}},\ldots,{\mathbf{1}}_{\\{g(x^{m})\geq
r^{m}\\}}\mathrel{\mathop{\mathchar
58\relax}}g\in\mathcal{G}\\}=\\{0,1\\}^{m},$

and no such $m$-tuple of points exists for $m>D$. This $D$ is called the VC-
subgraph dimension (or pseudo-dimension) of $\mathcal{G}$, and is denoted by
$\mathrm{vc}(\mathcal{G})$.

######  Lemma 4.4.

Let $\mathcal{G}$ be VC subgraph class of real-valued measurable functions
$g\mathrel{\mathop{\mathchar 58\relax}}{\mathbb{R}}^{d}\to[0,B]$. Let
$\boldsymbol{x}=(x^{1},\ldots,x^{N})$ be an arbitrary $N$-tuple of points in
${\mathbb{R}}^{d}$, and define the Rademacher average

$\displaystyle R_{N}(\mathcal{G};\boldsymbol{x})\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{N}{\mathbf{E}}\left[\sup_{g\in\mathcal{G}}\left|\sum^{N}_{i=1}\varepsilon^{i}g(x^{i})\right|\right],$

where $\varepsilon^{1},\ldots,\varepsilon^{N}$ are i.i.d. random variables
with $\mathbf{P}[\varepsilon^{i}=\pm 1]=\frac{1}{2}$. Then, for any
$N\geq\mathrm{vc}(\mathcal{G})$,

$\displaystyle R_{N}(\mathcal{G};\boldsymbol{x})\leq
cB\sqrt{\frac{\mathrm{vc}(\mathcal{G})\log N}{N}}$

for some universal constant $c$.

###  4.2 Proof of Theorem 3.1

Fix any $u\in\mathcal{U}$ and any $(\Sigma,\xi)$. Then

$\displaystyle\begin{split}&\phantom{{}={}}\|\mathsf{F}u-\mathsf{F}_{\Sigma,\xi}u\|_{\infty}\\\
&\qquad\leq\|\mathsf{F}u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F})u\|_{\infty}+\|\mathsf{F}_{\Sigma,\xi}u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F}_{\Sigma,\xi}\circ\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u\|_{\infty}\\\
&\qquad\qquad+\|(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F}_{\Sigma,\xi}\circ\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F})u\|_{\infty}\\\
&\qquad=\mathrel{\mathop{\mathchar 58\relax}}T_{1}+T_{2}+T_{3}.\end{split}$
(4.7)

Since $(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1})u=B_{k}u$, we can estimate
$T_{1}$ using Lemma 4.1 and Assumption 2.2:

$\displaystyle T_{1}$
$\displaystyle\leq\sup_{u\in\mathcal{U}}\|\mathsf{F}u-B_{k}(\mathsf{F}u,\cdot)\|_{\infty}=\sup_{y\in\mathcal{Y}}\|y-B_{k}(y,\cdot)\|_{\infty}\leq
2\omega_{\mathcal{Y}}\left(\frac{T}{\sqrt{k}}\right).$ (4.8)

For $T_{2}$, Lemma 4.2 gives

$\displaystyle T_{2}\leq
2\|\mathsf{F}_{\Sigma,\xi}\|_{\mathrm{Lip}}\omega_{u}\left(\frac{2T}{\sqrt{k}}\right)+2\omega_{\mathsf{F}_{\Sigma,\xi}B_{k-1}(u,\cdot)}\left(\frac{T}{\sqrt{k}}\right)$
(4.9)

Now, the system (3.1) has the form (4.2) with $f(x,u)=\sigma^{(n)}(Ax+Bu)$ and
$h=c$; thus, applying Lemma 4.3 to the i/o map $\mathsf{F}_{\Sigma,\xi}$ with
$\Sigma=(A,b,c)$, we get

$\displaystyle\|\mathsf{F}_{\Sigma,\xi}\|_{\mathrm{Lip}}\leq|c||b|e^{\|A\|T}\qquad\text{and}\qquad\omega_{\mathsf{F}_{\Sigma,\xi}B_{k-1}(u,\cdot)}(\delta)\leq\sqrt{n}|c|e^{\|A\|T}\delta.$

Substituting these estimates into (4.9) and invoking Assumption 2.1, we obtain

$\displaystyle T_{2}$ $\displaystyle\leq
2|c||b|e^{\|A\|T}\omega_{\mathcal{U}}\left(\frac{2T}{\sqrt{k}}\right)+2|c|Te^{\|A\|T}\sqrt{\frac{n}{k}}.$
(4.10)

Finally, using the fact that
$(\mathsf{S}^{*}_{k+1}\circ\mathsf{F}_{\Sigma,\xi}\circ\mathsf{S}^{*}_{k}\circ\mathsf{S}_{k})u=(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi}\circ\mathsf{S}_{k})u$, we can write

$\displaystyle T_{3}$ $\displaystyle=\|(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi}\circ\mathsf{S}_{k})u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F})u\|_{\infty}.$
(4.11)

Using (4.8), (4.10), and (4.11) in (4.7) and taking expectation with respect
to $\mu$, we obtain (3.2).

###  4.3 Proof of Theorem 3.2

For each pair $(\Sigma,\xi)$, define the function
$g_{\Sigma,\xi}\mathrel{\mathop{\mathchar
58\relax}}{\mathbb{R}}^{k}\times{\mathbb{R}}^{k+1}\to{\mathbb{R}}_{+}$
according to

$\displaystyle g_{\Sigma,\xi}(v,z)\mathrel{\mathop{\mathchar
58\relax}}=\max_{1\leq j\leq k}\left|(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi})v(t_{j})-\mathsf{S}^{*}_{k+1}z(t_{j})\right|,$

where $t_{j}=jT/k$. Let $\bar{\mu}$ denote the joint probability law of
$\mathsf{S}_{k}u$ and $\mathsf{S}_{k+1}y$ when $u\sim\mu$ and $y=\mathsf{F}u$.
Then $\bar{\mu}$ is a Borel probability measure on
${\mathbb{R}}^{k}\times{\mathbb{R}}^{k+1}$, and we can define the expected
risk

$\displaystyle\bar{\mathcal{L}}(\Sigma,\xi)\mathrel{\mathop{\mathchar
58\relax}}={\mathbf{E}}_{\bar{\mu}}[g_{\Sigma,\xi}(v,z)].$

Given the i/o data
$(u^{i},y^{i})\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mu$, the points
$(v^{i},z^{i})$ with $v^{i}=\mathsf{S}_{k}(u^{i})$ and
$z^{i}=\mathsf{S}_{k+1}(y^{i})$ are i.i.d. samples from $\bar{\mu}$, and our
learning procedure selects any minimizer $(\hat{\Sigma},\hat{\xi})$ of the
empirical risk

$\displaystyle\bar{\mathcal{L}}_{N}(\Sigma,\xi)\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{N}\sum^{N}_{i=1}g_{\Sigma,\xi}(v^{i},z^{i})$

among all $\Sigma=(A,b,c)$ and $\xi$ satisfying the constraint
$\|A\|,|b|,|c|,|\xi|\leq M$.

We next show that, with high probability, the excess risk
$\bar{\mathcal{L}}(\hat{\Sigma},\hat{\xi})-\bar{\mathcal{L}}^{*}$ is small,
where the minimum risk $\bar{\mathcal{L}}^{*}$ is defined in (3.5). By Lemma
4.3 applied to any i/o map $\mathsf{F}_{\Sigma,\xi}\in\mathcal{F}(M)$, and for
any $v\in C([0,T])$,

$\displaystyle\|B_{k}\mathsf{F}_{\Sigma,\xi}v\|_{\infty}\leq\|\mathsf{F}_{\Sigma,\xi}v\|_{\infty}\leq|c|\left(|\xi|+\sqrt{n}T\right)\leq
M(M+\sqrt{n}T).$

Moreover, for any $u\in\mathcal{U}$,

$\displaystyle\|(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F})u\|_{\infty}$
$\displaystyle=\|B_{k}(\mathsf{F}u,\cdot)\|_{\infty}\leq\|\mathsf{F}u\|_{\infty}\leq\gamma_{\mathsf{F}}(R).$

Therefore, we may assume without loss of generality that $0\leq
g_{\Sigma,\xi}(\cdot)\leq
M(M+\sqrt{n}T)+\gamma_{\mathsf{F}}(R)=\mathrel{\mathop{\mathchar 58\relax}}B$.
Then the usual ERM analysis (see, e.g., Cor. 6.1 in Hajek and Raginsky (2019))
guarantees that, with probability at least $1-\delta$,

$\displaystyle\bar{\mathcal{L}}(\hat{\Sigma},\hat{\xi})\leq\bar{\mathcal{L}}^{*}+4{\mathbf{E}}R_{N}(\mathcal{G})+B\sqrt{\frac{2\log(\frac{1}{\delta})}{N}},$

where
$R_{N}(\mathcal{G})=R_{N}(\mathcal{G};((v^{1},z^{1}),\ldots,(v^{N},z^{N})))$
is the Rademacher average of the function class
$\mathcal{G}\mathrel{\mathop{\mathchar
58\relax}}=\left\\{g_{\Sigma,\xi}\mathrel{\mathop{\mathchar
58\relax}}(\Sigma,\xi)\in\mathcal{F}(M)\right\\}$. Using Lemma 4.4, we then
see that

$\displaystyle\bar{\mathcal{L}}(\hat{\Sigma},\hat{\xi})\leq\bar{\mathcal{L}}^{*}+cB\sqrt{\frac{\mathrm{vc}(\mathcal{G})\log
N+\log(\frac{1}{\delta})}{N}}$ (4.12)

with probability at least $1-\delta$, provided
$N\geq\mathrm{vc}(\mathcal{G})$, where $c>0$ is an absolute constant and
$\mathrm{vc}(\mathcal{G})$ is the VC-subgraph dimension (or pseudo-dimension)
of $\mathcal{G}$.

#### Pseudo-dimension estimate.

For any $v\in{\mathbb{R}}^{k}$ and $z=(z_{0},\ldots,z_{k})^{\hbox{{\tiny
T}}}\in{\mathbb{R}}^{k+1}$, we can write

$\displaystyle g_{\Sigma,\xi}(v,z)=\max_{j\in[k]}|h^{(j)}_{\Sigma,\xi}(v,z)|,$

where

$\displaystyle h^{(j)}_{\Sigma,\xi}(v,z)\mathrel{\mathop{\mathchar
58\relax}}=\sum^{k}_{\ell=0}\big{(}y^{(\ell)}_{k,\Sigma,\xi}(v)-z_{\ell}\big{)}\frac{t^{\ell}_{j}}{\ell!},$

with $y^{(\ell)}_{k,\Sigma,\xi}(v)$ denoting the $\ell$th coordinate of
$Y_{k,\Sigma,\xi}(v)$. Let $\mathcal{H}^{(j)}$ denote the class of all
functions of the form $h^{(j)}_{\Sigma,\xi}$. We make the following two
claims:

  1. 1.

$\mathcal{H}^{(1)},\ldots,\mathcal{H}^{(k)}$ are VC-subgraph classes with the
same pseudo-dimension $d\leq 3n^{6}+5n^{3}\log_{2}k$;

  2. 2.

$\mathcal{G}$ is a VC-subgraph class with $\mathrm{vc}(\mathcal{G})\leq 2kd$.

We prove the second claim first. For each $j\in[k]$, consider the class
$\mathcal{C}^{(j)}$ of subsets of
${\mathbb{R}}^{k}\times{\mathbb{R}}^{k+1}\times{\mathbb{R}}$ of the form

$\displaystyle\left\\{(v,z,r)\in{\mathbb{R}}^{k}\times{\mathbb{R}}^{k+1}\times{\mathbb{R}}\mathrel{\mathop{\mathchar
58\relax}}|h^{(j)}_{\Sigma,\xi}(v,z)|\geq r\right\\}.$ (4.13)

Then evidently each set $\\{(v,z,r)\mathrel{\mathop{\mathchar
58\relax}}g_{\Sigma,\xi}(v,z)\geq r\\}$ is of the form $C^{(1)}\cup\ldots\cup
C^{(k)}$ with $C^{(j)}\in\mathcal{C}^{(j)}$. By the standard VC dimension
estimates (Giné and Nickl, 2016, Prop. 3.6.7), then,
$\mathrm{vc}(\mathcal{G})\leq\mathrm{vc}(\mathcal{H}^{(1)})+\ldots+\mathrm{vc}(\mathcal{H}^{(k)})$.
On the other hand, each set in (4.13) is itself the union of the sets
$\\{(v,z,r)\mathrel{\mathop{\mathchar 58\relax}}h^{(j)}_{\Sigma,\xi}(v,z)\geq
r\\}$ and $\\{(v,z,r)\mathrel{\mathop{\mathchar
58\relax}}h^{(j)}_{\Sigma,\xi}\leq-r\\}$, and therefore
$\mathrm{vc}(\mathcal{H}^{(j)})\leq 2d$. It remains to prove the first claim.

To that end, fix some $\Sigma=(A,b,c)$ and $\xi$ and let $\theta$ be a vector
of dimension $n^{2}+3n$ obtained by listing the entries of $A,b,c,\xi$ in some
fixed order. Then an induction argument together with the fact that the
function $\sigma(r)=\tanh r$ satisfies the identity
$\sigma^{\prime}(r)=1-\sigma^{2}(r)$ can be used to show that each function
$h^{(j)}_{\Sigma,\xi}$ can be written in the form

$\displaystyle
P(\sigma(R_{1}(\theta,v,z)),\ldots,\sigma(R_{n}(\theta,v,z)),\theta,v,z),$

where $P$ is a polynomial of degree at most $3k-1$ and each $R_{i}$ is a
polynomial of degree at most $2$ (Sontag, 1998). This implies, in turn, that

$\displaystyle\mathrm{vc}(\mathcal{H}^{(1)})=\ldots=\mathrm{vc}(\mathcal{H}^{(k)})\leq
3n^{6}+5n^{3}\log_{2}k$

(Sontag, 1998, Cor. 7). Substituting the resulting estimate of
$\mathrm{vc}(\mathcal{G})$ into (4.12), we see that

$\displaystyle\bar{\mathcal{L}}(\hat{\Sigma},\hat{\xi})\leq\bar{\mathcal{L}}^{*}+c(M(M+\sqrt{n}T)+\gamma_{\mathsf{F}}(R))\sqrt{\frac{k(n^{6}+n^{3}\log_{2}k)\log
N+\log(\frac{1}{\delta})}{N}}$ (4.14)

with probability at least $1-\delta$.

#### The final risk bound.

For any $\Sigma,\xi$ and any $u\in\mathcal{U}$,

$\displaystyle\|(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi}\circ\mathsf{S}_{k})u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F})u\|_{\infty}$
$\displaystyle=\|B_{k}(\mathsf{F}_{\Sigma,\xi}B_{k-1}u,\cdot)-B_{k}(\mathsf{F}u,\cdot)\|_{\infty}$
$\displaystyle=\sup_{t\in[0,T]}|B_{k}(\mathsf{F}_{\Sigma,\xi}B_{k-1}u,t)-B_{k}(\mathsf{F}u,t)|$
$\displaystyle\begin{split}&\leq\max_{1\leq j\leq
k}|B_{k}(\mathsf{F}_{\Sigma,\xi}B_{k-1}u,t_{j})-B_{k}(\mathsf{F}u,t_{j})|\\\
&\qquad+\sup_{t\in[0,T]}\min_{1\leq j\leq
k}|B_{k}(\mathsf{F}_{\Sigma,\xi}B_{k-1}u,t_{j})-B_{k}(\mathsf{F}_{\Sigma,\xi}B_{k-1}u,t)|\\\
&\qquad+\sup_{t\in[0,T]}\min_{1\leq j\leq
k}|B_{k}(\mathsf{F}u,t_{j})-B_{k}(\mathsf{F}u,t)|,\end{split}$

where the last two terms can be upper-bounded using the moduli of continuity
of $B_{k}(\mathsf{F}_{\Sigma,\xi}B_{k-1}u,\cdot)$ and
$B_{k}(\mathsf{F}u,\cdot)$, which, by Lemma 4.1, are upper-bounded by twice
the moduli of continuity of $\mathsf{F}_{\Sigma,\xi}B_{k-1}u$ and
$\mathsf{F}u$, respectively. Thus, using Assumption 2.2 and Lemma 4.3, we get

$\displaystyle\begin{split}&\|(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi}\circ\mathsf{S}_{k})u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F})u\|_{\infty}\\\
&\qquad\leq\max_{1\leq j\leq
k}|B_{k}(\mathsf{F}_{\Sigma,\xi}B_{k-1}u,t_{j})-B_{k}(\mathsf{F}u,t_{j})|+\frac{2}{k}\sqrt{n}MTe^{MT}+2\omega_{\mathcal{Y}}\left(\frac{T}{k}\right).\end{split}$

Taking expectation of both sides with respect to $\mu$ yields, for any
$(\Sigma,\xi)\in\mathcal{F}(M)$,

$\displaystyle{\mathbf{E}}_{\mu}\left[\|(\mathsf{S}^{*}_{k+1}\circ
Y_{k,\Sigma,\xi}\circ\mathsf{S}_{k})u-(\mathsf{S}^{*}_{k+1}\circ\mathsf{S}_{k+1}\circ\mathsf{F})u\|_{\infty}\right]\leq\bar{\mathcal{L}}(\Sigma,\xi)+\frac{2}{k}\sqrt{n}MTe^{MT}+2\omega_{\mathcal{Y}}\left(\frac{T}{k}\right).$

Using this and (4.14) in the bound of Theorem 3.1, we get the statement of the
theorem.

## Acknowledgments

The work of J. Hanson and M. Raginsky was supported in part by the Illinois
Institute for Data Science and Dynamical Systems (iDS2), an NSF HDR TRIPODS
institute, under award CCF-1934986 and by the Center for Advanced Electronics
through Machine Learning (CAEML) NSF I/UCRC award CNS-16-24811.

## References

  * Chua (1980) Leon O. Chua.  Device modeling via basic nonlinear circuit elements.  _IEEE Transactions on Circuits and Systems_ , CAS-27(11):1014–1044, November 1980.
  * DeVore and Lorentz (1993) Ronald A. DeVore and George G. Lorentz.  _Constructive Approximation_.  Springer-Verlag, New York, 1993.
  * Farrell et al. (2020) Max H. Farrell, Tengyuan Liang, and Sanjog Misra.  Deep neural networks for estimation and inference.  _Econometrica_ , 2020.  URL https://arxiv.org/abs/1809.09953.  To appear.
  * Funahashi and Nakamura (1993) Ken-Ichi Funahashi and Yuichi Nakamura.  Approximation of dynamical systems by continuous time recurrent neural networks.  _Neural Networks_ , 6(6):801 – 806, 1993.
  * Giné and Nickl (2016) Evarist Giné and Richard Nickl.  _Mathematical Foundations of Infinite-Dimensional Statistical Models_.  Cambridge University Press, 2016.
  * Hajek and Raginsky (2019) Bruce Hajek and Maxim Raginsky.  ECE 543: Statistical Learning Theory.  University of Illinois lecture notes, 2019.  URL http://maxim.ece.illinois.edu/teaching/SLT/SLT.pdf.
  * Hanson and Raginsky (2020) Joshua Hanson and Maxim Raginsky.  Universal simulation of stable dynamical systems by recurrent neural nets.  In _Proc. 2nd Conference on Learning for Dynamics and Control_ , volume 120 of _Proceedings of Machine Learning Research_ , pages 384–392, 2020.
  * Hirsch and Smale (1974) Morris W. Hirsch and Stephen Smale.  Nonautonomous equations and differentiability of flows.  In _Differential equations, dynamical systems, and linear algebra_ , chapter 15, pages 296–303. Academic Press, 1974.
  * Li (2000) Zhongkai Li.  Bernstein polynomials and modulus of continuity.  _Journal of Approximation Theory_ , 102:171–174, 2000.
  * Sontag (1992) Eduardo D. Sontag.  Neural nets as systems models and controllers.  In _Proc. Seventh Yale Workshop on Adaptive and Learning Systems_ , pages 73–79, 1992.
  * Sontag (1998) Eduardo D. Sontag.  A learning result for continuous-time recurrent neural networks.  _Systems and Control Letters_ , 34:151–158, 1998.