Framing RNN as a kernel method: A neural ODE approach
Year: 2022

Abstract:
Building on the interpretation of a recurrent neural network (RNN) as a continuoustime neural differential equation, we show, under appropriate conditions, that the solution of a RNN can be viewed as a linear function of a specific feature set of the input sequence, known as the signature. This connection allows us to frame a RNN as a kernel method in a suitable reproducing kernel Hilbert space. As a consequence, we obtain theoretical guarantees on generalization and stability for a large class of recurrent networks. Our results are illustrated on simulated datasets.

heading: 1 Introduction
text: Recurrent neural networks (RNN) are among the most successful methods for modeling sequential data. They have achieved state-of-the-art results in difficult problems such as natural language processing (e.g., Mikolov et al., 2010; Collobert et al., 2011) or speech recognition (e.g., Hinton et al., 2012; Graves et al., 2013). This class of neural networks has a natural interpretation in terms of (discretization of) ordinary differential equations (ODE), which casts them in the field of neural ODE (Chen et al., 2018). This observation has led to the development of continuous-depth models for handling irregularly-sampled time-series data, including the ODE-RNN model (Rubanova et al., 2019), GRU-ODE-Bayes (De Brouwer et al., 2019), or neural CDE models (Kidger et al., 2020; Morrill et al., 2020a). In addition, the time-continuous interpretation of RNN allows to leverage the rich theory of differential equations to develop new recurrent architectures (Chang et al., 2019; Herrera et al., 2020; Erichson et al., 2021), which are better at learning long-term dependencies.
On the other hand, the development of kernel methods for deep learning offers theoretical insights on the functions learned by the networks (Cho and Saul, 2009; Belkin et al., 2018; Jacot et al., 2018). Here, the general principle consists in defining a reproducing kernel Hilbert space (RKHS)—that is, a function class H —, which is rich enough to describe the architectures of networks. A good example is the construction of Bietti and Mairal (2017, 2019), who exhibit an RKHS for convolutional neural networks. This kernel perspective has several advantages. First, by separating the representation of the data from the learning process, it allows to study invariances of the representations learned by the network. Next, by reducing the learning problem to a linear one in H , generalization bounds can be more easily obtained. Finally, the Hilbert structure of H provides a natural metric on neural networks, which can be used for example for regularization (Bietti et al., 2019).
∗Equal contribution
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Contributions. By taking advantage of the neural ODE paradigm for RNN, we show that RNN are, in the continuous-time limit, linear predictors over a specific space associated with the signature of the input sequence (Levin et al., 2013). The signature transform, first defined by Chen (1958) and central in rough path theory (Lyons et al., 2007; Friz and Victoir, 2010), summarizes sequential inputs by a graded feature set of their iterated integrals. Its natural environment is a tensor space that can be endowed with an RKHS structure (Király and Oberhauser, 2019). We exhibit general conditions under which classical recurrent architectures such as feedforward RNN, Gated Recurrent Units (GRU, Cho et al., 2014), or Long Short-Term Memory networks (LSTM, Hochreiter and Schmidhuber, 1997), can be framed as a kernel method in this RKHS. This enables us to provide generalization bounds for RNN as well as stability guarantees via regularization. The theory is illustrated with some experimental results.
Related works. The neural ODE paradigm was first formulated by Chen et al. (2018) for residual neural networks. It was then extended to RNN in several articles, with a focus on handling irregularly sampled data (Rubanova et al., 2019; Kidger et al., 2020) and learning long-term dependencies (Chang et al., 2019). The signature transform has recently received the attention of the machine learning community (Levin et al., 2013; Kidger et al., 2019; Liao et al., 2019; Toth and Oberhauser, 2020; Fermanian, 2021) and, combined with deep neural networks, has achieved state-of-the-art performance for several applications (Yang et al., 2016, 2017; Perez Arribas, 2018; Wang et al., 2019; Morrill et al., 2020b). Király and Oberhauser (2019) use the signature transform to define kernels for sequential data and develop fast computational methods. The connection between continuous-time RNN and signatures has been pointed out by Lim (2021) for a specific model of stochastic RNN. Deriving generalization bounds for RNN is an active research area (Zhang et al., 2018; Akpinar et al., 2019; Tu et al., 2019). By leveraging the theory of differential equations, our approach encompasses a large class of RNN models, ranging from feedforward RNN to LSTM. This is in contrast with most existing generalization bounds, which are architecture-dependent. Close to our point of view is the work of Bietti and Mairal (2017) for convolutional neural networks.
Mathematical context. We place ourselves in a supervised learning setting. The input data is a sample of n i.i.d. vector-valued sequences {x(1), . . . ,x(n)}, where x(i) = (x(i)1 , . . . , x (i) T ) ∈ (Rd)T , T ≥ 1. The outputs of the learning problem can be either labels (classification setting) or sequences (sequence-to-sequence setting). Even if we only observe discrete sequences, each x(i) is mathematically considered as a regular discretization of a continuous-time process X(i) ∈ BV c([0, 1],Rd), where BV c([0, 1],Rd) is the space of continuous functions from [0, 1] to Rd of finite total variation. Informally, the total variation of a process corresponds to its length. Formally, for any [s, t] ⊂ [0, 1], the total variation of a process X ∈ BV c([0, 1],Rd) on [s, t] is defined by
‖X‖TV ;[s,t] = sup (t0,...,tk)∈Ds,t k∑ j=1 ‖Xtj −Xtj−1‖,
where Ds,t denotes the set of all finite partitions of [s, t] and ‖ · ‖ the Euclidean norm. We therefore have that x(i)j = X (i) j/T , 1 ≤ j ≤ T , where X (i) t := X
(i)(t). We make two assumptions on the processes X(i). First, they all begin at zero, and second, their lengths are bounded by L ∈ (0, 1). These assumptions are not too restrictive, since they amount to data translation and normalization, common in practice. Accordingly, we denote by X the subset of BV c([0, 1],Rd) defined by
X = { X ∈ BV c([0, 1],Rd) |X0 = 0 and ‖X‖TV ;[0,1] ≤ L } and assume therefore that X(1), . . . , X(n) are i.i.d. according to some X ∈ X . The norm on all spaces Rm, m ≥ 1, is always the Euclidean one. Observe that assuming that X ∈X implies that, for any t ∈ [0, 1], ‖Xt‖ = ‖Xt −X0‖ ≤ ‖X‖TV ;[0,1] ≤ L.
Recurrent neural networks. Classical RNN are defined by a sequence of hidden states h1, . . . , hT ∈ Re, where, for x = (x1, . . . , xT ) a generic data sample,
h0 = 0 and hj+1 = f(hj , xj+1) for 0 ≤ j ≤ T − 1.
At each time step 1 ≤ j ≤ T , the output of the network is zj = ψ(hj), where ψ is a linear function. In the present article, we rather consider the following residual version, which is a natural adaptation
of classical RNN in the neural ODE framework (see, e.g., Yue et al., 2018):
h0 = 0 and hj+1 = hj + 1
T f(hj , xj+1) for 0 ≤ j ≤ T − 1. (1)
The simplest choice for the function f is the feedforward model, say fRNN, defined by
fRNN(h, x) = σ(Uh+ V x+ b), (2)
where σ is an activation function, U ∈ Re×e and V ∈ Re×d are weight matrices, and b ∈ Re is the bias. The function fRNN, equipped with a smooth activation σ (such as the logistic or hyperbolic tangent functions), will be our leading example throughout the paper. However, the GRU and LSTM models can also be rewritten under the form (1), as shown in Appendix A.1. Thus, model (1) is flexible enough to encompass most recurrent networks used in practice.
Overview. Section 2 is devoted to framing RNN as linear functions in a suitable RKHS. We start by embedding iteration (1) into a continuous-time model, which takes the form of a controlled differential equation (CDE). This allows, after introducing the signature transform, to define the appropriate RKHS, and, in turn, to show that model (1) boils down, in the continuous-time limit, to a linear problem on the signature. This framework is used in Section 3 to derive generalization bounds and stability guarantees. We provide some experiments in Section 4 before discussing our results in Section 5. All proofs are postponed to the supplementary material.

heading: 2 Framing RNN as a kernel method
text: Roadmap. First, we quantify the difference between the discrete recurrent network (1) and its continuous-time counterpart (Proposition 1). Then, we rewrite the corresponding ODE as a CDE (Proposition 2). Under appropriate conditions, Proposition 4 shows that the solution of this equation is a linear function of the signature of the driving process. Importantly, these assumptions are valid for a feedforward RNN, as stated by Proposition 5. We conclude in Theorem 1.

heading: 2.1 From discrete to continuous time
text: Recall that h0, . . . , hT denote the hidden states of the RNN (1), and let H : [0, 1] → Re be the solution of the ODE dHt = f(Ht, Xt)dt, H0 = h0. (3) By bounding the difference between Hj/T and hj , the following proposition shows how to pass from discrete to continuous time, provided f satisfies the following assumption:
(A1) The function f is Lipschitz continuous in h and x, with Lipschitz constants Kh and Kx. We let Kf = max(Kh,Kx).
Proposition 1. Assume that (A1) is verified. Then there exists a unique solution H to (3) and, for any 0 ≤ j ≤ T ,
‖Hj/T − hj‖ ≤ c1 T ,
where c1 = KfeKf ( L + sup
‖h‖≤M,‖x‖≤L ‖f(h, x)‖eKf
) and M = sup
‖x‖≤L ‖f(h0, x)‖eKf . Moreover,
for any t ∈ [0, 1], ‖Ht‖ ≤M .
Then, following Kidger et al. (2020), we show that the ODE (3) can be rewritten under the form of a CDE. At the cost of increasing the dimension of the hidden state from e to e+ d, this allows us to reframe model (3) as a linear model in dX , in the sense that X has been moved ‘outside’ of f . Proposition 2. Assume that (A1) is verified. Let H : [0, 1] → Re be the solution of (3), and let X̄ : [0, 1]→ Rd+1 be the time-augmented process X̄t = (X>t , 1−L2 t)
>. Then there exists a tensor field F : Rē → Rē×d̄, ē = e+ d, d̄ = d+ 1, such that if H̄ : [0, 1]→ Rē is the solution of the CDE
dH̄t = F(H̄t)dX̄t, H̄0 = (H > 0 , X > 0 ) >, (4)
then its first e coordinates are equal to H .
Equation (4) can be better understood by the following equivalent integral equation:
H̄t = H̄0 + ∫ t 0 F(H̄u)dX̄u,
where the integral should be understood as Riemann-Stieljes integral (Friz and Victoir, 2010, Section I.2). Thus, the output of the RNN can be approximated by the solution of the CDE (4), and, according to Proposition 1, the approximation error is O(1/T). Example 1. Consider fRNN as in (2). If σ is Lipschitz continuous with constant Kσ, then, for any h1, h2 ∈ Re, x1, x2 ∈ Rd,
‖fRNN(h1, x1)− fRNN(h2, x1)‖ = ‖σ(Uh1 + V x1 + b)− σ(Uh2 + V x1 + b)‖ ≤ Kσ‖U‖op‖h1 − h2‖,
where ‖ · ‖op denotes the operator norm—see Appendix A.3. Similarly, ‖f(h1, x1)− f(h1, x2)‖ ≤ Kσ‖V ‖op‖x1 − x2‖. Thus, assumption (A1) is satisfied. The tensor field FRNN of Proposition 2 corresponding to this network is defined for any h̄ ∈ Rē by
FRNN(h̄) = ( 0e×d 21−Lσ(Wh̄+ b) Id×d 0d×1 ) , where W = (U V ) ∈ Re×ē. (5)

heading: 2.2 The signature
text: An essential ingredient towards our construction is the signature of a continuous-time process, which we briefly present here. We refer to Chevyrev and Kormilitzin (2016) for a gentle introduction and to Lyons et al. (2007); Levin et al. (2013) for details.
Tensor Hilbert spaces. We denote by (Rd)⊗k the kth tensor power of Rd with itself, which is a Hilbert space of dimension dk. The key space to define the signature and, in turn, our RKHS, consists in infinite square-summable sequences of tensors of increasing order:
T = { a = (a0, . . . , ak, . . . ) ∣∣∣ ak ∈ (Rd)⊗k, ∞∑ k=0 ‖ak‖2(Rd)⊗k <∞ } . (6)
Endowed with the scalar product 〈a, b〉T := ∑∞ k=0〈ak, bk〉(Rd)⊗k , T is a Hilbert space, as shown in Appendix A.4. Definition 1. Let X ∈ BV c([0, 1],Rd). For any t ∈ [0, 1], the signature of X on [0, t] is defined by S[0,t](X) = (1,X1[0,t], . . . ,X k [0,t], . . . ), where, for each k ≥ 1,
Xk[0,t] = k! ∫ · · · ∫
0≤u1<···<uk≤t
dXu1 ⊗ · · · ⊗ dXuk ∈ (Rd)⊗k.
Although this definition is technical, the signature should simply be thought of as a feature map that embeds a bounded variation process into an infinite-dimensional tensor space. The signature has several good properties that make it a relevant tool for machine learning (e.g., Levin et al., 2013; Chevyrev and Kormilitzin, 2016; Fermanian, 2021). In particular, under certain assumptions, S(X) characterizes X up to translations and reparameterizations, and has good approximation properties. We also highlight that fast libraries exist for computing the signature (Reizenstein and Graham, 2020; Kidger and Lyons, 2021).
The expert reader is warned that this definition differs from the usual one by the normalization of Xk[0,t] by k!, which is more adapted to our context. In the sequel, for any index (i1, . . . , ik) ⊂ {1, . . . , d}k, S
(i1,...,ik) [0,t] (X) denotes the term associated with the coordinates (i1, . . . , ik) of X k [0,t]. When the signature is taken on the whole interval [0, 1], we simply write S(X), S(i1,...,ik)(X), and Xk. Example 2. Let X be the d-dimensional linear path defined by Xt = (a1 + b1t, . . . , ad + bdt)>, ai, bi ∈ R. Then S(i1,...,ik)(X) = bi1 · · · bik and Xk = b⊗k.
The next proposition, which ensures that S[0,t](X̄) ∈ T , is an important step. Proposition 3. Let X ∈ X and X̄t = (X>t , 1−L2 t)
> as in Proposition 2. Then, for any t ∈ [0, 1], ‖S[0,t](X̄)‖T ≤ 2(1− L)−1.
The signature kernel. By taking advantage of the structure of Hilbert space of T , it is natural to introduce the following kernel:
K : X ×X → R (X,Y ) 7→ 〈S(X̄), S(Ȳ )〉T ,
which is well defined according to Proposition 3. We refer to Király and Oberhauser (2019) for a general presentation of kernel methods with signatures and to Cass et al. (2020) for a kernel trick. The RKHS associated with K is the space of functions
H = { ξα : X → R | ξα(X) = 〈α, S(X̄)〉T , α ∈ T } , (7)
with scalar product 〈ξα, ξβ〉H = 〈α, β〉T (see, e.g., Schölkopf and Smola, 2002).

heading: 2.3 From the CDE to the signature kernel
text: An important property of signatures is that the solution of the CDE (4) can be written, under certain assumptions, as a linear function of the signature of the driving process X . This operation can be thought of as a Taylor expansion for CDE. More precisely, let us rewrite (4) as
dHt = F(Ht)dXt = d∑ i=1 F i(Ht)dX i t , (8)
where Xt = (X1t , . . . , X d t ) >, F : Re → Re×d, and F i : Re → Re are the columns of F—to avoid heavy notation, we momentarily write e, d, H , and X instead of ē, d̄, H̄ , and X̄ . Throughout, the bold notation is used to distinguish tensor fields and vector fields. We recall that a vector field F : Re → Re or a tensor field F : Re → Re×d are said to be smooth if each of their coordinates is C∞. Definition 2. Let F,G : Re → Re be smooth vector fields and denote by J(·) the Jacobian matrix. Their differential product is the smooth vector field F ? G : Re → Re defined, for any h ∈ Re, by
(F ? G)(h) = e∑ j=1 ∂G ∂hj (h)Fj(h) = J(G)(h)F (h).
In differential geometry, F ?G is simply denoted by FG. Since the ? operation is not associative, we take the convention that it is evaluated from right to left, i.e., F 1 ? F 2 ? F 3 := F 1 ? (F 2 ? F 3).
Taylor expansion. Let H be the solution of (8), where F is assumed to be smooth. We now show that H can be written as a linear function of the signature of X , which is the crucial step to embed the RNN in the RKHS H . The step-N Taylor expansion of H (Friz and Victoir, 2008) is defined by
HNt = H0 + N∑ k=1 1 k! ∑ 1≤i1,...,ik≤d S (i1,...,ik) [0,t] (X)F i1 ? · · · ? F ik(H0).
Throughout, we let Λk(F) = sup
‖h‖≤M,1≤i1,...,ik≤d ‖F i1 ? · · · ? F ik(h)‖.
Example 3. Let F = FRNN defined by (5) with an identity activation. Then, for any h̄ ∈ Rē, 1 ≤ i ≤ d+ 1, F iRNN(h̄) = Wih̄+ bi, where bi is the (i+ d)th vector of the canonical basis of Rē, and
Wi = 0ē×ē, Wd+1 =
( 2
1−LW 0d×ē
) , and bd+1 = ( 2
1−Lb 0d
) .
The vector fields F iRNN are then affine, J(F i RNN) = Wi, and the iterated star products have a simple expression: for any 1 ≤ i1, . . . , ik ≤ d, F i1RNN ? · · · ? F ik RNN(h̄) = Wik · · ·Wi2(Wi1 h̄+ bi1).
The next proposition shows that the step-N Taylor expansion HN is a good approximation of H . Proposition 4. Assume that the tensor field F is smooth. Then, for any t ∈ [0, 1], N ≥ 1,
‖Ht −HNt ‖ ≤ dN+1
(N + 1)! ΛN+1(F). (9)
Thus, provided that ΛN (F) is not too large, the right-hand side of (9) converges to zero, hence
Ht = H0 + ∞∑ k=1 1 k! ∑ 1≤i1,...,ik≤d S (i1,...,ik) [0,t] (X)F i1 ? · · · ? F ik(H0). (10)
We conclude from the above representation that the solution H of (8) is in fact a linear function of the signature of X . A natural concern is to know whether the upper bound of Proposition 4 vanishes with N for standard architectures. This property is encapsulated in the following more general assumption:
(A2) The tensor field F is smooth and ∞∑ k=0 (dk k! Λk(F) )2 <∞.
Clearly, if (A2) is verified, then the right-hand side of (9) converges to 0. The next proposition states formally the conditions under which (A2) is verified for FRNN. It is further illustrated in Figure 1, which shows that the convergence is fast with two common activation functions. We let ‖σ‖∞ = sup‖h‖≤M,‖x‖≤L ‖σ(Uh+V x+b)‖ and ‖σ(k)‖∞ = sup‖h‖≤M,‖x‖≤L ‖σ(k)(Uh+V x+b)‖, where σ(k) is the derivative of order k of σ. Proposition 5. Let FRNN be defined by (5). If σ is the identity function, then (A2) is satisfied. In the general case, (A2) holds if σ is smooth and there exists a > 0 such that, for any k ≥ 0,
‖σ(k)‖∞ ≤ ak+1k! and ‖W‖F < 1− L 8a2d , (11)
where ‖ · ‖F is the Frobenius norm. Moreover, for any N ≥ 1,
ΛN (FRNN) ≤ √ 2a (8a2‖W‖F
1− L
)N−1 N ! .
The proof of Proposition 5, based on the manipulation of higher-order derivatives of tensor fields, is highly non-trivial. We highlight that the conditions on σ are mild and verified for common smooth activations. For example, they are verified for the logistic function (with a = 2) and for the hyperbolic tangent function (with a = 4)—see Appendix A.5. The second inequality of (11) puts a constraint on the norm of the weights, and can be regarded as a radius of convergence for the Taylor expansion.
Putting everything together. We now have all the elements at hand to embed the RNN into the RKHS H . To fix the idea, we assume in this paragraph that we are in a ±1 classification setting. In other words, given an input sequence x, we are interested in the final output zT = ψ(hT ) ∈ R, where hT is the solution of (1). The predicted class is 2 · 1(zT > 0)− 1. By Propositions 1 and 2, zT is approximated by the first e coordinates of the solution of the CDE (4), which outputs a Re+d-valued process H̄ . According to Proposition 4, H̄ is a linear function of the signature of the time-augmented process X̄ . Thus, on top of H̄ , it remains to successively apply the projection Proj on the e first coordinates followed by the linear function ψ to obtain an element of the RKHS H . This mechanism is summarized in the following theorem. Theorem 1. Assume that (A1) and (A2) are verified. Then there exists a function ξα ∈H such that
|zT − ξα(X)| ≤ ‖ψ‖op c1 T , (12)
where ξα(X) = 〈α, S(X̄)〉T and X̄t = (X>t , 1−L2 t) >. We have α = (αk)∞k=0, where each αk ∈ (Rd)⊗k is defined by
α (i1,...,ik) k =
1 k! ψ ◦ Proj
( F i1 ? · · · ? F ik(H̄0) ) .
Moreover, ‖α‖2T ≤ ‖ψ‖2op ∑∞ k=0 ( dk k! Λk(F) )2 .
We conclude that in the continuous-time limit, the output of the network can be interpreted as a scalar product between the signature of the (time-augmented) process X̄ and an element of T . This interpretation is important for at least two reasons: (i) it facilitates the analysis of generalization of RNN by leveraging the theory of kernel methods, and (ii) it provides new insights on regularization strategies to make RNN more robust. These points will be explored in the next section. Finally, we stress that the approach works for a large class of RNN, such as GRU and LSTM. The derivation of conditions (A1) and (A2) beyond the feedforward RNN is left for future work.

heading: 3 Generalization and regularization

heading: 3.1 Generalization bounds
text: Learning procedure. A first consequence of framing a RNN as a kernel method is that it gives natural generalization bounds under mild assumptions. In the learning setup, we are given an i.i.d. sample Dn of n random pairs of observations (x(i),y(i)) ∈ (Rd)T × Y , where x(i) = (x
(i) 1 , . . . , x (i) T ). We distinguish the binary classification problem, where Y = {−1, 1}, from the sequential prediction problem, where Y = (Rp)T and y(i) = (y(i)1 , . . . , y (i) T ). The RNN is assumed to be parameterized by θ ∈ Θ ⊂ Rq, where Θ is a compact set. To clarify the notation, we use a θ subscript whenever a quantity depends on θ (e.g., fθ for f , etc.). In line with Section 2, it is assumed that the tensor field Fθ associated with fθ satisfies (A1) and (A2), keeping in mind that Proposition 5 guarantees that these requirements are fulfilled by a feedforward recurrent network with a smooth activation function.
Let gθ : (Rd)T → Y denote the output of the recurrent network. The parameter θ is fitted by empirical risk minimization using a loss function ` : Y × Y → R+. The theoretical and empirical risks are respectively defined, for any θ ∈ Θ, by
R(θ) = E[`(y, gθ(x))] and R̂n(θ) = 1
n n∑ i=1 ` ( y(i), gθ(x (i)) ) ,
where the expectation E is evaluated with respect to the distribution of the generic random pair (x,y). We let θ̂n ∈ argmin θ∈ΘR̂n(θ) and aim at upper bounding P ( ygθ̂n(x) ≤ 0|Dn ) in the classification regime (Theorem 2) and R(θ̂n) in the sequential regime (Theorem 3). To reach this goal, our strategy is to approximate the RNN by its continuous version and then use the RKHS machinery of Section 2.
Binary classification. In this context, the network outputs a real number gθ(x) = ψ(hT ) ∈ R and the predicted class is 2 · 1(gθ(x) > 0)− 1. The loss ` : {−1, 1} × R → R+ is assumed to satisfy the assumptions of Bartlett and Mendelson (2002, Theorem 7), that is, for any y ∈ {−1, 1}, z ∈ R, `(y, z) = φ(yz), where φ(u) ≥ 1(u ≤ 0), and φ is Lipschitz-continuous with constant K`. For example, the logistic loss satisfies such assumptions. We let ξαθ ∈H be the function of Theorem 1 that approximates the RNN with parameter θ. Thus, zT ≈ ξαθ (X̄) = 〈αθ, S(X̄)〉T , up to a O(1/T) term. Theorem 2. Assume that for all θ ∈ Θ, (A1) and (A2) are verified. Assume, in addition, that there exists a constant B > 0 such that for any θ ∈ Θ, ‖ξαθ‖H ≤ B. Then with probability at least 1− δ,
P ( ygθ̂n(x) ≤ 0|Dn ) ≤ R̂n(θ̂n) + c2 T + 8BK` (1− L) √ n + 2BK` 1− L
√ log(1/δ)
2n , (13) where c2 = K` supθ ( ‖ψ‖opKfθeKfθ ( L+ ‖fθ‖∞eKfθ )) .
Close to our result are the bounds obtained by Zhang et al. (2018), Tu et al. (2019), and Chen et al. (2020). The main difference is that the term in 1/T does not usually appear, since it comes from the Euler discretization error, whereas the speed in 1/√n is the same. For instance, Chen et al. (2020) show that, under some assumptions, the excess risk is of order √ de+ e2TαK`n
−1/2. We refer to Section 5 for further discussion on the dependency of the different bounds to the parameter T . The take-home message is that the detour by continuous-time neural ODE provides a theoretical framework adapted to RNN, at the modest price of an additional O(1/T) term. Moreover, we note that the bound (13) is ‘simple’ and holds under mild conditions for a large class of RNN. More precisely, for any recurrent network of the form (1), provided (A1) and (A2) are satisfied, then (13) is valid with constants c2 and B depending on the architecture. Such constants are given below in the example of a feedforward RNN. We stress that Theorem 2 can be extended without significant effort to the multi-class classification task, with an appropriate choice of loss function. Example 4. Take a feedforward RNN with logistic activation, and Θ = {(W, b, ψ) | ‖W‖F ≤ KW < (1− L)/32d, ‖b‖ ≤ Kb, ‖ψ‖op ≤ Kψ}. Then, Proposition 5 states that (A2) is satisfied and, with Theorem 1, ensures that
sup θ∈Θ ‖ξαθ‖H ≤
√ 2Kψ(1− L)
1− L− 32dKW := B, Kfθ = max(‖U‖op, ‖V ‖op), and ‖fθ‖∞ = 1.
Sequence-to-sequence learning. We conclude by showing how to extend both the RKHS embedding of Theorem 1 and the generalization bound of Theorem 2 to the setting of sequence-to-sequence learning. In this case, the output of the network is a sequence
gθ(x) = (z1, . . . , zT ) ∈ (Rp)T .
An immediate extension of Theorem 1 ensures that there exist p elements α1,θ, . . . , αp,θ ∈ T such that, for any 1 ≤ j ≤ T ,∥∥zj − (〈α1,θ, S[0,j/T ](X̄)〉T , . . . , 〈αp,θ, S[0,j/T ](X̄)〉T )>∥∥ ≤ ‖ψ‖op c1
T . (14)
The properties of the signature guarantee that S[0,j/T ](X) = S(X̃[j]) where X̃[j] is the process equal to X̄ on [0, j/T ] and then constant on [j/T , 1]—see Appendix A.6. With this trick, we have, for any 1 ≤ ` ≤ p, 〈α`,θ, S[0,j/T ](X̄)〉T = 〈α`,θ, S(X̃[j])〉T , so that we are back in H . Observe that the only difference with (12) is that we consider vector-valued sequential outputs, which requires to introduce the process X̃[j], but that the rationale is exactly the same.
We let ` : (Rp)T × (Rp)T → R+ be the L2 distance, that is, for any y = (y1, . . . , yT ), y′ = (y′1, . . . , y ′ T ), `(y,y ′) = 1T ∑T j=1 ‖yj − y′j‖2. It is assumed that y takes its values in a compact subset of Rq , i.e., there exists Ky > 0 such that ‖yj‖ ≤ Ky . Theorem 3. Assume that for all θ ∈ Θ, (A1) and (A2) are verified. Assume, in addition, that there exists a constant B > 0 such that for any 1 ≤ ` ≤ p, θ ∈ Θ, ‖ξα`,θ‖H ≤ B . Then with probability at least 1− δ,
R(θ̂n) ≤ R̂n(θ̂n) + c3 T
+ 4pc4B(1− L)−1√
n +
√ 2c5 log(1/δ)
n , (15)
where c3 = sup θ
( c1,θ + ‖ψ‖op‖fθ‖∞ ) + 2 √ pB(1 − L)−1 + 2Ky, c4 = B(1 − L)−1 + Ky, and
c5 = 4pB(1− L)−1c4 +K2y .

heading: 3.2 Regularization and stability
text: In addition to providing a sound theoretical framework, framing deep learning in an RKHS provides a natural norm, which can be used for regularization, as shown for example in the context of convolutional neural networks by Bietti et al. (2019). This regularization ensures stability of predictions, which is crucial in particular in a small sample regime or in the presence of adversarial examples (Gao et al., 2018; Ko et al., 2019). In our binary classification setting, for any inputs x,x′ ∈ (Rd)T , by the Cauchy-Schwartz inequality, we have
‖zT − z′T ‖ ≤ 2‖ψ‖op‖ c1 T + ‖ξαθ (X̄)− ξαθ (X̄ ′)‖ ≤ 2‖ψ‖op‖ c1 T + ‖ξαθ‖H ‖S(X̄)− S(X̄ ′)‖T .
If x and x′ are close, so are their associated continuous processes X and X ′ (which can be approximated for example by taking a piecewise linear interpolation), and so are their signatures. The term ‖S(X̄) − S(X̄ ′)‖T is therefore small (Friz and Victoir, 2010, Proposition 7.66). Therefore, when T is large, we see that the magnitude of ‖ξαθ‖H determines how close the predictions are. A natural training strategy to ensure stable predictions, for the types of networks covered in the present article, is then to penalize the problem by minimizing the loss R̂n(θ) + λ‖ξαθ‖2H . From a computational point of view, it is possible to compute the norm in H , up to a truncation at N of the Taylor expansion, which we know by Proposition 4 to be reasonable. It remains that computing this norm is a non-trivial task, and implementing smart surrogates is an interesting problem for the future. Note however that computing the signature of the data is not necessary for this regularization strategy.

heading: 4 Numerical illustrations
text: This section is here for illustration purposes. Our objective is not to achieve competitive performance, but rather to illustrate the theoretical results. We refer to Appendix D for implementation details.
Convergence of the Taylor expansion towards the solution of the ODE. We illustrate Proposition 4 on a toy example. The process X is a 2-dimensional spiral, and we take feedforward RNN with 2 hidden units. Repeating this procedure with 103 uniform random weight initializations, we observe in Figure 1a that the signature approximation converges exponentially fast in N . As seen in Figure 1b, the rate of convergence depends in particular on the norm of the weight matrices, as predicted by Proposition 5. However, condition (11) seems to be over-restrictive, since convergence happens even for weights with norm larger than the bound (we have 1/(8a2d) ' 0.01 here).
Adversarial robustness. We illustrate the penalization proposed in Section 3.2 on a toy task that consists in classifying the rotation direction of 2-dimensional spirals. We take a feedforward RNN with 32 hidden units and hyperbolic tangent activation. It is trained on 50 examples, with and without penalization, for 200 epochs. Once trained, the RNN is tested on adversarial examples, generated with the projected gradient descent algorithm with Frobenius norm (Madry et al., 2018), which modifies test examples to maximize the error while staying in a ball of radius ε. We observe in Figure 2 that adding the penalization seems to make the network more stable.
Comparison of the trained networks. The evolution of the Frobenius norm of the weights ‖W‖F and the RKHS norm ‖ξαθ‖H during training is shown in Figure 3. This points out that the penalization, which forces the RNN to keep a small norm in H , leads indeed to learning different weights than the non-penalized RNN. The results also suggest that the Frobenius and RKHS norms are decoupled, since both networks have Frobenius norms of similar magnitude but very different RKHS norms. The figures show one random run, but we observe similar qualitative behavior on others.

heading: 5 Discussion and conclusion
text: Role of the discretization procedure. The starting point of the paper was motivated by the fact that the classical residual RNN formulation coincides with an Euler discretization of the ODE (3).
This choice of discretization translates into a 1/T term in Theorems 2 and 3. However, we could have considered higher-order discretization schemes, such as Runge-Kutta schemes, for which the discretization error decreases as 1/Tp. Such schemes correspond to alternative architectures, which were already proposed by Wang and Lin (1998), among others. At the limit, we could also consider directly the continuous model (3), as proposed by Chen et al. (2018), in which case the discretization error term vanishes. Of course, such an option requires to be able to sample the continuous-time data at arbitrary times.
Long-term stability. RNN are known to be poor at learning long-term dependencies (Bengio et al., 1993; Hochreiter and Schmidhuber, 1997). This is reflected in the literature by performance bounds increasing in T , which is not the case of our results (13) and (15), seemingly indicating that we fail to capture this phenomenon. This apparent paradox is related to our assumption that the total variation of X is bounded. Indeed, if a time series is observed for a long time, then its total variation may become large. In this case, it is no longer valid to assume that ‖X‖TV is bounded by L. In other words, in our context, the parameter encapsulating the notion of “long-term” is not T but the regularity of X measured by its total variation. Note that the choice of defining X on [0, 1] and not another interval [0, U ] is arbitrary and does not carry any meaning on the problem of learning long-term dependencies. A thorough analysis of these questions is an interesting research direction for future work.
Radius of convergence. The assumptions ‖X‖TV ;[0,1] ≤ L < 1 and ‖W‖F ≤ KW < (1−L)/32d can be seen as radii of convergence of the Taylor expansion (10). They allow using the Taylor approximation—which is of a local nature—to prove a global result, the RKHS embedding. In return, the condition on the Frobenius norm of the weights puts restrictions on the admissible parameters of the neural network. However, this bound can be improved, in particular by considering more exotic norms, which we did not explicit for clarity purposes.
Conclusion. By bringing together the theory of neural ODE, the signature transform, and kernel methods, we have shown that a recurrent network can be framed in the continuous-time limit as a linear function in a well-chosen RKHS. In addition to giving theoretical insights on the function learned by the network and providing generalization guarantees, this framing suggests regularization strategies to obtain more robust RNN. We have only scratched the surface of the potentialities of leveraging this theory to practical applications, which is a subject of its own and will be tackled in future work.

heading: Acknowledgements
text: Authors thank T. Lévy for his inputs on the Picard-Lindelöf theorem and N. Doumèche for fruitful discussion. A. Fermanian has been supported by a grant from Région Île-de-France and P. Marion by a stipend from Corps des Mines.

Review summaries:
The authors consider a continuous data/neural ode model for recurrent neural networks. They prove convergence of the discrete model to the continuum model as the length of input sequences increases (assuming that sequences are obtained by finer and finer discretizations of the same process). They then show that the neural ODE can be represented alternatively as a controlled ODE, and that a solution of the controlled ODE can be expressed in terms of the signature of the input process. This allows the authors to equate the RNN to a kernel method for a specific kernel. Finally, they use Rademacher complexity bounds for the RKHS to estimate the generalization properties of RNNs in terms of the length of the input sequence/fineness of discretization. They use the continuum model to design a penalty term which improves stability as seen in numerical experiments.

This paper presents an approach to interpret RNN (in the continuous-time setting) as a kernel method. The main idea is to write the RNN ODE as a controlled differential equation, and through a signature transform of the (time-augmented) input process it is shown that the RNN output can be viewed as a linear function of the signature. Hence, RNN viewed this way is now amenable to classical analysis from Kernel methods. In particular, a generalization bound is proved in which the complexity measure is based on the defined RKHS norm of the target. Moreover, the approach inspires some regularization strategies for RNNs and the paper explores some of this via toy examples.

This paper provides a rigorous framework to connect learning with RNNs and kernel methods, exploiting the idea that the output of a RNN can be viewed as a linear function of the signature of the input sequence in the continuous-time setting. Working within the framework (embedding the RNN into a RKHS), generalization bounds and stability guarantees are provided. Experiments are also provided to support the theory.