figurec

# The Recurrent Neural Tangent Kernel

Sina Alemohammad, Zichao Wang, Randall Balestriero, Richard G. Baraniuk
Department of Electrical and Computer Engineering
Rice University
<EMAIL_ADDRESS>
###### Abstract

The study of deep neural networks (DNNs) in the infinite-width limit, via the
so-called neural tangent kernel (NTK) approach, has provided new insights into
the dynamics of learning, generalization, and the impact of initialization.
One key DNN architecture remains to be kernelized, namely, the recurrent
neural network (RNN). In this paper we introduce and study the Recurrent
Neural Tangent Kernel (RNTK), which provides new insights into the behavior of
overparametrized RNNs. A key property of the RNTK should greatly benefit
practitioners is its ability to compare inputs of different length. To this
end, we characterize how the RNTK weights different time steps to form its
output under different initialization parameters and nonlinearity choices. A
synthetic and 56 real-world data experiments demonstrate that the RNTK offers
significant performance gains over other kernels, including standard NTKs,
across a wide array of data sets.

##  1 Introduction

The overparameterization of modern deep neural networks (DNNs) has resulted in
not only remarkably good generalization performance on unseen data [32, 31, 7]
but also guarantees that gradient descent learning can find the global minimum
of their highly nonconvex loss functions [15, 2, 1, 39, 5]. From these
successes, a natural question arises: What happens when we take
overparameterization to the limit by allowing the width of a DNN’s hidden
layers to go to infinity? Surprisingly, the analysis of such an (impractical)
DNN becomes analytically tractable. Indeed, recent work has shown that the
training dynamics of (infinite-width) DNNs under gradient flow is captured by
a constant kernel called the _Neural Tangent Kernel_ (NTK) that evolves
according to a linear ordinary differential equation (ODE) [24, 29, 4].

Every DNN architecture and parameter initialization produces a distinct NTK.
The original NTK was derived from the Multilayer Perceptron (MLP)[24] and was
soon followed by kernels derived from Convolutional Neural Networks (CNTK) [4,
35], Residual DNNs [23], and Graph Convolutional Neural Networks (GNTK) [14].
In [37], a general strategy to obtain the NTK of any architecture is provided.

In this paper, we extend the NTK concept to the important class of
overparametrized Recurrent Neural Networks (RNNs), a fundamental DNN
architecture for processing sequential data. We show that RNN in its infinite-
width limit converges to a kernel that we dub the Recurrent Neural Tangent
Kernel (RNTK). The RNTK provides high performance for various machine learning
tasks, and an analysis of the properties of the kernel provides useful
insights into the behavior of RNNs in the following overparametrized regime.
In particular, we derive and study the RNTK to answer the following
theoretical questions:

Q: Can the RNTK extract long-term dependencies between two data sequences?
RNNs are known to underperform at learning long-term dependencies due to the
gradient vanishing or exploding [8]. Attempted ameliorations have included
orthogonal weights [3, 25, 20] and gating such as in Long Short-Term Memory
(LSTM) [21] and Gated Recurrent Unit (GRU) [11] RNNs. We demonstrate that the
RNTK can detect long-term dependencies with proper initialization of the
hyperparameters, and moreover, we show how the dependencies are extracted
through time via different hyperparameter choices.

Q: Do the recurrent weights of the RNTK reduce its representation power
compared to other NTKs? An attractive property of an RNN that is shared by the
RNTK is that it can deal with sequences of different lengths via weight
sharing through time. This enables the reduction of the number of learnable
parameters and thus more stable training at the cost of reduced representation
power. We prove the surprising fact that employing tied vs. untied weights in
an RNN does not impact the analytical form of the RNTK.

Q: Does the RNTK generalize well? A recent study has revealed that the use of
an SVM classifier with the NTK, CNTK, and GNTK kernels outperforms other
classical kernel-based classifiers and trained finite DNNs on small data sets
(typically fewer than 5000 training samples) [27, 4, 6, 14]. We extend these
results to RNTKs to demonstrate that the RNTK outperforms a variety of classic
kernels, NTKs and finite RNNs for time series data sets in both classification
and regression tasks. Carefully designed experiments with data of varying
lengths demonstrate that the RNTK’s performance accelerates beyond other
techniques as the difference in lengths increases. Those results extend the
empirical observations from [4, 6, 14, 27] into finite DNNs, NTK, CNTK, and
GNTK comparisons by observing that their performance-wise ranking depends on
the employed DNN architecture.

We summarize our contributions as follows:

[C1] We derive the analytical form for the RNTK of an overparametrized RNN at
initialization using rectified linear unit (ReLU) and error function (erf)
nonlinearities for arbitrary data lengths and number of layers (Section 3.1).

[C2] We prove that the RNTK remains constant during (overparametrized) RNN
training and that the dynamics of training are simplified to a set of ordinary
differential equations (ODEs) (Section 3.2).

[C3] When the input data sequences are of equal length, we show that the RNTKs
of weight-tied and weight-untied RNNs converge to the same RNTK (Section 3.3).

[C4] Leveraging our analytical formulation of the RNTK, we empirically
demonstrate how correlations between data at different times are weighted by
the function learned by an RNN for different sets of hyperparameters. We also
offer practical suggestions for choosing the RNN hyperparameters for deep
information propagation through time (Section 3.4).

[C5] We demonstrate that the RNTK is eminently practical by showing its
superiority over classical kernels, NTKs, and finite RNNs in exhaustive
experiments on time-series classification and regression with both synthetic
and 56 real-world data sets (Section 4).

##  2 Background and Related Work

Notation. We denote $[n]=\\{1,\dots,n\\}$, and ${\bm{I}}_{d}$ the identity
matrix of size $d$. $[{\bm{A}}]_{i,j}$ represents the $(i,j)$-th entry of a
matrix, and similarly $[{\bm{a}}]_{i}$ represents the $i$-th entry of a
vector. We use $\phi(\cdot):\mathbb{R}\rightarrow\mathbb{R}$ to represent the
activation function that acts coordinate wise on a vector and $\phi^{\prime}$
to denote its derivative. We will often use the rectified linear unit (ReLU)
$\phi(x)=\mathrm{max}(0,x)$ and error function (erf)
$\phi(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-z^{2}}dz$ activation functions.
$\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$ represents the multidimensional
Gaussian distribution with the mean vector ${\bm{\mu}}$ and the covariance
matrix ${\bm{\Sigma}}$.

Recurrent Neural Networks (RNNs). Given an input sequence data
${\bm{x}}=\\{{\bm{x}}_{t}\\}_{t=1}^{T}$ of length $T$ with data at time $t$,
${\bm{x}}_{t}\in\mathbb{R}^{m}$, a simple RNN [17] performs the following
recursive computation at each layer $\ell$ and each time step $t$

$\displaystyle{\bm{g}}^{\left(\ell,t\right)}({\bm{x}})={\bm{W}}^{\left(\ell\right)}{\bm{h}}^{\left(\ell,t-1\right)}({\bm{x}})+{\bm{U}}^{\left(\ell\right)}{\bm{h}}^{\left(\ell-1,t\right)}({\bm{x}})+{\bm{b}}^{\left(\ell\right)},\hskip
20.0pt{\bm{h}}^{\left(\ell,t\right)}({\bm{x}})$
$\displaystyle=\phi\left({\bm{g}}^{\left(\ell,t\right)}({\bm{x}})\right),$

where ${\bm{W}}^{(\ell)}\in\mathbb{R}^{n\times n}$,
${\bm{b}}^{(\ell)}\in\mathbb{R}^{n}$ for $\ell\in[L]$,
${\bm{U}}^{(1)}\in\mathbb{R}^{n\times m}$ and
${\bm{U}}^{(\ell)}\in\mathbb{R}^{n\times n}$ for $\ell\geq 2$ are the RNN
parameters. ${\bm{g}}^{\left(\ell,t\right)}({\bm{x}})$ is the pre-activation
vector at layer $\ell$ and time step $t$, and
${\bm{h}}^{\left(\ell,t\right)}({\bm{x}})$ is the after-activation (hidden
state). For the input layer $\ell=0$, we define
${\bm{h}}^{\left(0,t\right)}({\bm{x}}):={\bm{x}}_{t}$.
${\bm{h}}^{(\ell,0)}({\bm{x}})$ as the initial hidden state at layer $\ell$
that must be initialized to start the RNN recursive computation.

The output of an $L$-hidden layer RNN with linear read out layer is achieved
via

$\displaystyle f_{\theta}({\bm{x}})$
$\displaystyle={\bm{V}}{\bm{h}}^{(L,T)}({\bm{x}}),$

where ${\bm{V}}\in\mathbb{R}^{d\times n}$. Figure 1 visualizes an RNN unrolled
through time.

${\bm{h}}^{(2,2)}$${\bm{h}}^{(2,1)}({\bm{x}})$${\bm{h}}^{(2,3)}({\bm{x}})$${\bm{h}}^{(1,2)}({\bm{x}})$${\bm{h}}^{(1,1)}({\bm{x}})$${\bm{h}}^{(1,3)}({\bm{x}})$${\bm{W}}^{(2)}$${\bm{W}}^{(2)}$${\bm{W}}^{(1)}$${\bm{W}}^{(1)}$${\bm{W}}^{(1)}$${\bm{U}}^{(2)}$${\bm{U}}^{(2)}$${\bm{U}}^{(2)}$${\bm{h}}^{(2,0)}({\bm{x}})$${\bm{h}}^{(1,0)}({\bm{x}})$${\bm{W}}^{(2)}$${\bm{W}}^{(1)}$${\bm{x}}_{1}$${\bm{x}}_{2}$${\bm{x}}_{3}$${\bm{U}}^{(1)}$${\bm{U}}^{(1)}$${\bm{U}}^{(1)}$

Figure 1:  Visualization of a simple RNN that highlights a cell (purple), a
layer (red) and the initial hidden state of each layer (green). (Best viewed
in color.)

Neural Tangent Kernel (NTK). Let $f_{\theta}({\bm{x}})\in\mathbb{R}^{d}$ be
the output of a DNN with parameters $\theta$. For two input data sequences
${\bm{x}}$ and ${\bm{x}}^{\prime}$, the NTK is defined as [24]

$\displaystyle\widehat{\Theta}_{s}({\bm{x}},{\bm{x}}^{\prime})=\langle\nabla_{\theta_{s}}f_{\theta_{s}}({\bm{x}}),\nabla_{\theta_{s}}f_{\theta_{s}}({\bm{x}}^{\prime})\rangle,$

where $f_{\theta_{s}}$ and $\theta_{s}$ are the network output and parameters
during training at time s.111We use $s$ to denote time here, since $t$ is used
to index the time steps of the RNN inputs.  Let $\mathcal{X}$ and
$\mathcal{Y}$ be the set of training inputs and targets,
$\ell(\widehat{y},y):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{+}$
be the loss function, and
$\mathcal{L}=\frac{1}{{|\mathcal{X}|}}\sum_{({\bm{x}},{\bm{y}})\in\mathcal{X}\times\mathcal{Y}}\ell(f_{\theta_{s}}({\bm{x}}),{\bm{y}})$
be the the empirical loss. The evolution of the parameters $\theta_{s}$ and
output of the network $f_{\theta_{s}}$ on a test input using gradient descent
with infinitesimal step size (a.k.a gradient flow) with learning rate $\eta$
is given by

$\displaystyle\frac{\partial\theta_{s}}{\partial s}$
$\displaystyle=-\eta\nabla_{\theta_{s}}f_{\theta_{s}}(\mathcal{X})^{T}\nabla_{f_{\theta_{s}}(\mathcal{X})}\mathcal{L}$
(1) $\displaystyle\frac{\partial f_{\theta_{s}}({\bm{x}})}{\partial s}$
$\displaystyle=-\eta\nabla_{\theta_{s}}f_{\theta_{s}}({\bm{x}})\nabla_{\theta_{s}}f_{\theta_{s}}(\mathcal{X})^{T}\nabla_{f_{\theta_{s}}(\mathcal{X})}\mathcal{L}=-\eta\widehat{\Theta}_{s}({\bm{x}},\mathcal{X})\nabla_{f_{\theta_{s}}(\mathcal{X})}\mathcal{L}.$
(2)

Generally, $\widehat{\Theta}_{s}({\bm{x}},{\bm{x}}^{\prime})$, hereafter
referred to as the empirical NTK, changes over time during training, making
the analysis of the training dynamics difficult. When $f_{\theta_{s}}$
corresponds to an infinite-width MLP, [24] showed that
$\widehat{\Theta}_{s}({\bm{x}},{\bm{x}}^{\prime})$ converges to a limiting
kernel at initialization and stays constant during training, i.e.,

$\displaystyle\underset{n\rightarrow\infty}{\mathrm{lim}}\widehat{\Theta}_{s}({\bm{x}},{\bm{x}}^{\prime})=\underset{n\rightarrow\infty}{\mathrm{lim}}\widehat{\Theta}_{0}({\bm{x}},{\bm{x}}^{\prime}):=\Theta({\bm{x}},{\bm{x}}^{\prime})\,\,\,\,\forall
s\,,$

which is equivalent to replacing the outputs of the DNN by their first-order
Taylor expansion in the parameter space [29]. With a mean-square error (MSE)
loss function, the training dynamics in (1) and (2) simplify to a set of
linear ODEs, which coincides with the training dynamics of kernel ridge
regression with respect to the NTK when the ridge term goes to zero. A nonzero
ridge regularization can be conjured up by adding a regularization term
$\frac{\lambda^{2}}{2}\|\theta_{s}-\theta_{0}\|_{2}^{2}$ to the empirical loss
[22].

##  3 The Recurrent Neural Tangent Kernel

We are now ready to derive the RNTK. We first prove the convergence of an RNN
at initialization to the RNTK in the infinite-width limit and discuss various
insights it provides. We then derive the convergence of an RNN after training
to the RNTK. Finally, we analyze the effects of various hyperparameter choices
on the RNTK. Proofs of all of our results are provided in the Appendices.

###  3.1 RNTK for an Infinite-Width RNN at Initialization

First we specify the following parameter initialization scheme that follows
previous work on NTKs [24], which is crucial to our convergence results:

$\displaystyle{\bm{W}}^{(\ell)}=\frac{\sigma_{w}^{\ell}}{\sqrt{n}}{\mathbf{W}}^{(\ell)},\,\,\,\,{\bm{U}}^{(1)}=\frac{\sigma_{u}^{1}}{\sqrt{m}}{\mathbf{U}}^{(1)},\,\,\,\,{\bm{U}}^{(\ell)}=\frac{\sigma_{u}^{\ell}}{\sqrt{n}}{\mathbf{U}}^{(\ell)}(\ell\geq
2),\,\,\,\,{\bm{V}}=\frac{\sigma_{v}}{\sqrt{n}}{\mathbf{V}},\,\,\,\,{\bm{b}}^{(\ell)}=\sigma_{b}{\mathbf{b}}^{(\ell)}\,,$
(3)

where

$\displaystyle[{\mathbf{W}}^{\ell}]_{i,j},\,\,[{\mathbf{U}}^{(\ell)}]_{i,j},\,\,[{\mathbf{V}}]_{i,j},\,\,[{\mathbf{b}}^{(\ell)}]_{i}\sim\mathcal{N}(0,1)\,.$
(4)

We will refer to (3) and (4) as the NTK initialization. The choices of the
hyperparameters $\sigma_{w}$, $\sigma_{u}$, $\sigma_{v}$ and $\sigma_{b}$ can
significantly impact RNN performance, and we discuss them in detail in Section
3.4. For the initial (at time $t=0$) hidden state at each layer $\ell$, we set
${\bm{h}}^{(\ell,0)}({\bm{x}})$ to an i.i.d. copy of
$\mathcal{N}(0,\sigma_{h})$ [34] . For convenience, we collect all of the
learnable parameters of the RNN into
$\theta=\mathrm{vect}\big{[}\\{\\{{\mathbf{W}}^{(\ell)},{\mathbf{U}}^{(\ell)},{\mathbf{b}}^{(\ell)}\\}_{\ell=1}^{L},{\mathbf{V}}\\}\big{]}$.

The derivation of the RNTK at initialization is based on the correspondence
between Gaussian initialized, infinite-width DNNs and Gaussian Processes
(GPs), known as the DNN-GP. In this setting every coordinate of the DNN output
tends to a GP as the number of units/neurons in the hidden layer (its width)
goes to infinity. The corresponding DNN-GP kernel is computed as

$\displaystyle\mathcal{K}({\bm{x}},{\bm{x}}^{\prime})=\mathop{{\mathbb{E}}}_{\theta\sim\mathcal{N}}\big{[}[f_{\theta}({\bm{x}})]_{i}\bm{\cdot}[f_{\theta}({\bm{x}}^{\prime})]_{i}\big{]},\hskip
2.84526pt\forall i\in[d].$ (5)

First introduced for a single-layer, fully-connected neural network by [30],
recent works on NTKs have extended the results for various DNN architectures
[28, 16, 33, 19, 36], where in addition to the output, all pre-activation
layers of the DNN tends to a GPs in the infinite-width limit. In the case of
RNNs, each coordinate of the RNN pre-activation
${\bm{g}}^{(\ell,t)}({\bm{x}})$ converges to a centered GP depending on the
inputs with kernel

$\displaystyle\Sigma^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})=\mathop{{\mathbb{E}}}_{\theta\sim\mathcal{N}}\big{[}[{\bm{g}}^{(\ell,t)}({\bm{x}})]_{i}\bm{\cdot}[{\bm{g}}^{(\ell,t^{\prime})}{({\bm{x}}^{\prime})}]_{i}\big{]}\,\,\,\forall
i\in[n].$ (6)

As per [35], the gradients of random infinite-width DNNs computed during
backpropagation are also Gaussian distributed. In the case of RNNs, every
coordinate of the vector
${\bm{\delta}}^{(\ell,t)}({\bm{x}}):=\sqrt{n}\big{(}\nabla_{{\bm{g}}^{(\ell,t)}({\bm{x}})}f_{\theta}({\bm{x}})\big{)}$
converges to a GP with kernel

$\displaystyle\Pi^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})=\mathop{{\mathbb{E}}}_{\theta\sim\mathcal{N}}\big{[}[{\bm{\delta}}^{(\ell,t)}({\bm{x}})]_{i}\bm{\cdot}[{\bm{\delta}}^{(\ell,t^{\prime})}({\bm{x}}^{\prime})]_{i}\big{]}\,\,\,\forall
i\in[n].$ (7)

Both convergences occur independently of the coordinate index $i$ and for
inputs of possibly different lengths, i.e., $T\neq T^{\prime}$. With (6) and
(7), we now prove that an infinite-width RNN at initialization converges to
the limiting RNTK.

###### Theorem 1

Let ${\bm{x}}$ and ${\bm{x}}^{\prime}$ be two data sequences of potentially
different lengths $T$ and $T^{\prime}$, respectively. Without loss of
generality, assume that $T\leq T^{\prime}$, and let $\tau:=T^{\prime}-T$. Let
$n$ be the number of units in the hidden layers, the empirical RNTK for an
$L$-layer RNN with NTK initialization converges to the following limiting
kernel as $n\rightarrow\infty$

$\displaystyle\underset{n\rightarrow\infty}{\mathrm{lim}}\widehat{\Theta}_{0}({\bm{x}},{\bm{x}}^{\prime})=\Theta({\bm{x}},{\bm{x}}^{\prime})=\Theta^{(L,T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\otimes{\bm{I}}_{d}\,,$
(8)

where

$\displaystyle\Theta^{(L,T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})=\left(\sum_{\ell=1}^{L}\sum_{t=1}^{T}\left(\Pi^{(\ell,t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Sigma^{(\ell,t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})\right)\right)+\mathcal{K}({\bm{x}},{\bm{x}}^{\prime})\,,$
(9)

with $\mathcal{K}({\bm{x}},{\bm{x}}^{\prime})$,
$\Sigma^{(\ell,t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$, and
$\Pi^{(\ell,t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$ defined in (5)–(7).

Remarks. Theorem 1 holds generally for any two data sequences, including
different lengths ones. This highlights the RNTK’s ability to produce a
similarity measure $\Theta({\bm{x}},{\bm{x}}^{\prime})$ even if the inputs are
of different lengths, without resorting to heuristics such as zero padding the
inputs to the to the max length of both sequences. Dealing with data of
different length is in sharp contrast to common kernels such as the classical
radial basis functions, polynomial kernels, and current NTKs. We showcase this
capability below in Section 4.

To visualize Theorem 1, we plot in the left plot in Figure 2 the convergence
of a single layer, sufficiently wide RNN to its RNTK with the two simple
inputs ${\bm{x}}=\\{1,-1,1\\}$ of length 3 and ${\bm{x}}^{\prime}=\\{{\rm
cos}(\alpha),{\rm sin}(\alpha)\\}$ of length 2, where $\alpha=[0,2\pi]$. For
an RNN with a sufficiently large hidden state ($n=1000$), we see clearly that
it converges to the RNTK ($n=\infty$).

$\,\,\,\,\,\,\,\,\,\,\widehat{\Theta}_{0}(x,x^{\prime})$

$\alpha$

$\,\,\,\,\,\,\,\,\,\,\,\,\,\,\log(\|\frac{\Theta-\widehat{\Theta}_{0}}{\Theta}\|)$

$\log(n)$

Figure 2:  Empirical demonstration of a wide, single-layer RNN converging to
its limiting RNTK. Left: convergence for a pair of different-length inputs
${\bm{x}}=\\{1,-1,1\\}$ and
${\bm{x}}^{\prime}=\\{\cos(\alpha),\sin(\alpha)\\}$, with varying
$\alpha=[0,2\pi]$. The vertical axis corresponds to the RNTK values for
different values of $\alpha$. Right: convergence of weight-tied and weight-
untied single layer RNN to the same limiting RNTK with increasing width
(horizontal axis). The vertical axis corresponds to the average of the log-
normalized error between the empirical RNTK computed using finite RNNs and the
RNTK for 50 Gaussian normal signals of length $T=5$.

RNTK Example for a Single-Layer RNN. We present a concrete example of Theorem
1 by showing how to recursively compute the RNTK for a single-layer RNN; thus
we drop the layer index for notational simplicity. We compute and display the
RNTK for the general case of a multi-layer RNN in Appendix B.3. To compute the
RNTK $\Theta^{(T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$, we need to compute
the GP kernels $\Sigma^{(t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$ and
$\Pi^{(t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$. We first define the operator
$\mathrm{V}_{\phi}\big{[}{\bm{K}}\big{]}$ that depends on the nonlinearity
$\phi(\cdot)$ and a positive semi-definite matrix
${\bm{K}}\in\mathbb{R}^{2\times 2}$

$\displaystyle\mathrm{V}_{\phi}\big{[}{\bm{K}}\big{]}={\mathbb{E}}[\phi({\textnormal{z}}_{1})\bm{\cdot}\phi({\textnormal{z}}_{2})],\qquad({\textnormal{z}}_{1},{\textnormal{z}}_{2})\sim\mathcal{N}(0,{\bm{K}})\,.$
(10)

Following [35], we obtain the analytical recursive formula for the GP kernel
$\Sigma^{(t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$ for a single layer RNN as

$\displaystyle\Sigma^{(1,1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\sigma_{h}^{2}1_{({\bm{x}}={\bm{x}}^{\prime})}+\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{1},{\bm{x}}^{\prime}_{1}\rangle+\sigma_{b}^{2}$
(11) $\displaystyle\Sigma^{(t,1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{1}\rangle+\sigma_{b}^{2}$
$\displaystyle t>1$ (12)
$\displaystyle\Sigma^{(1,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{1},{\bm{x}}^{\prime}_{t^{\prime}}\rangle+\sigma_{b}^{2}$
$\displaystyle t^{\prime}>1$ (13)
$\displaystyle\Sigma^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\big{]}+\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{t^{\prime}}\rangle+\sigma_{b}^{2}\hskip
2.84526pt$ $\displaystyle t,t^{\prime}>1$ (14)
$\displaystyle\mathcal{K}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(T+1,T^{\prime}+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]},$
(15)

where

$\displaystyle{\bm{K}}^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\left[{\begin{array}[]{cc}\Sigma^{(t-1,t-1)}({\bm{x}},{\bm{x}})&\Sigma^{(t-1,t^{\prime}-1)}({\bm{x}},{\bm{x}}^{\prime})\\\
\Sigma^{(t-1,t^{\prime}-1)}({\bm{x}},{\bm{x}}^{\prime})&\Sigma^{(t^{\prime}-1,t^{\prime}-1)}({\bm{x}}^{\prime},{\bm{x}}^{\prime})\\\
\end{array}}\right].$ (18)

Similarly, we obtain the analytical recursive formula for the GP kernel
$\Pi^{(t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$ as

$\displaystyle\Pi^{(T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(T+1,T+\tau+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}$
(19) $\displaystyle\Pi^{(t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(t+1,t+\tau+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}\Pi^{(t+1,t+1+\tau)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle t\in[T-1]$ (20)
$\displaystyle\Pi^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=0$ $\displaystyle t^{\prime}-t\neq\tau.$ (21)

For $\phi=\mathrm{ReLU}$ and $\phi=\mathrm{erf}$, we provide analytical
expressions for $\mathrm{V}_{\phi}\big{[}{\bm{K}}\big{]}$ and
$\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}\big{]}$ in Appendix B.5. These
yield an explicit formula for the RNTK that enables fast and point-wise kernel
evaluations. For other activation functions, one can apply the Monte Carlo
approximation to obtain $\mathrm{V}_{\phi}\big{[}{\bm{K}}\big{]}$ and
$\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}\big{]}$ [33].

###  3.2 RNTK for an Infinite-Width RNN during Training

We prove that an infinitely-wide RNN, not only at initialization but also
during gradient descent training, converges to the limiting RNTK at
initialization.

###### Theorem 2

Let $n$ be the number of units of each RNN’s layer. Assume that
$\Theta(\mathcal{X},\mathcal{X})$ is positive definite on $\mathcal{X}$ such
that $\lambda_{\mathrm{min}}(\Theta(\mathcal{X},\mathcal{X}))>0$. Let
$\eta^{*}:=2\big{(}\lambda_{\mathrm{min}}(\Theta(\mathcal{X},\mathcal{X}))+\lambda_{\mathrm{max}}(\Theta(\mathcal{X},\mathcal{X})\big{)}\big{)}^{-1}$.
For an $L$-layer RNN with NTK initialization as in (3), (4) trained under
gradient flow (recall (1) and (2)) with $\eta<\eta^{*}$, we have with high
probability

$\displaystyle\underset{s}{\mathrm{sup}}\frac{\|\theta_{s}-\theta_{0}\|_{2}}{\sqrt{n}},\underset{s}{\mathrm{sup}}\|\widehat{\Theta}_{s}(\mathcal{X},\mathcal{X})-\widehat{\Theta}_{0}(\mathcal{X},\mathcal{X})\|_{2}=\mathcal{O}\left(\frac{1}{\sqrt{n}}\right).$

Remarks. Theorem 2 states that the training dynamics of an RNN in the
infinite-width limit as in (1), (2) are governed by the RNTK derived from the
RNN at its initialization. Intuitively, this is due to the NTK initialization
(3), (4) which positions the parameters near a local minima, thus minimizing
the amount of update that needs to be applied to the weights to obtain the
final parameters.

$\frac{s(t)}{\max_{t}s(t)}$

$\frac{s(t)}{\max_{t}s(t)}$

$t$

$t$

$t$

$t$

Figure 3: Per time step $t$ (horizontal axis) sensitivity analysis (vertical
axis) of the RNTK for the ReLU (top row) and erf (bottom row) activation
functions for various weight noise hyperparameters. We also experiment with
different RNTK hyperparameters in each of the subplots, given by the subplot
internal legend. Clearly, the ReLU (top-row) provides a more stable kernel
across time steps (highlighted by the near constant sensitivity through time).
On the other hand, erf (bottom row) sees a more erratic behavior either
focusing entirely on early time-steps or on the latter ones.

###  3.3  RNTK for an Infinite-Width RNN Without Weight Sharing

We prove that, in the infinite-width limit, an RNN without weight sharing
(untied weights), i.e., using independent new weights
${\mathbf{W}}^{(\ell,t)}$, ${\mathbf{U}}^{(\ell,t)}$ and
${\mathbf{b}}^{(\ell,t)}$ at each time step $t$, converges to the same RNTK as
an RNN with weight sharing (tied weights). First, recall that it is a common
practice to use weight-tied RNNs, i.e., in layer $\ell$, the weights
${\mathbf{W}}^{(\ell)}$, ${\mathbf{U}}^{(\ell)}$ and ${\mathbf{b}}^{(\ell)}$
are the same across all time steps $t$. This practice conserves memory and
reduces the number of learnable parameters. We demonstrate that, when using
untied-weights, the RNTK formula remains unchanged.

###### Theorem 3

For inputs of the same length, an RNN with untied weights converges to the
same RNTK as an RNN with tied weights in the infinite-width
($n\rightarrow\infty$) regime.

Remarks. Theorem 3 implies that weight-tied and weight-untied RNNs have
similar behaviors in the infinite-width limit. It also suggests that existing
results on the simpler, weight-untied RNN setting may be applicable for the
more general, weight-tied RNN. The plot on the right side of Figure 2
empirically demonstrates the convergence of both the weight-tied and weight-
untied RNNs to the RNTK with increasing hidden layer size $n$; moreover, the
convergence rates are similar.

###  3.4 Insights into the Roles of the RNTK’s Hyperparameters

Our analytical form for the RNTK is fully determined by a small number of
hyperparameters, which contains the various weight variances collected into
$\mathcal{S}=\\{\sigma_{w},\sigma_{u},\sigma_{b},\sigma_{h}\\}$ and the
activation function.222From (11) to (20) we emphasize that $\sigma_{v}$ merely
scales the RNTK and does not change its overall behavior. In standard
supervised-learning settings, one often performs cross-validation to select
the hyperparameters. However, since kernel methods become computationally
intractable for large datasets, we seek a more computationally friendly
alternative to cross-validation. Here we conduct a novel exploratory analysis
that provides new insights into the impact of the RNTK hyperparameters on the
RNTK output and suggests a simple method to select them a priori in a
deliberate manner.

To visualize the role of the RNTK hyperparameters, we introduce the
sensitivity $s(t)$ of the RNTK of two input sequences ${\bm{x}}$ and
${\bm{x}}^{\prime}$ with respect to the input ${\bm{x}}_{t}$ at time $t$

$\displaystyle
s(t)=\|\nabla_{{\bm{x}}_{t}}\Theta({\bm{x}},{\bm{x}}^{\prime})\|_{2}\,.$ (22)

Here, $s(t)$ indicates how sensitive the RNTK is to the data at time $t$,
i.e., ${\bm{x}}_{t}$, in presence of another data sequence
${\bm{x}}^{\prime}$. Intuitively, large/small $s(t)$ indicates that the RNTK
is relatively sensitive/insensitive to the input ${\bm{x}}_{t}$ at time $t$.

The sensitivity is crucial to understanding to which extent the RNTK
prediction is impacted by the input at each time step. In the case where some
time indices have a small sensitivity, then any input variation in those
corresponding times will not alter the RNTK output and thus will produce a
metric that is invariant to those changes. This situation can be beneficial or
detrimental based on the task at hand. Ideally, and in the absence of prior
knowledge on the data, one should aim to have a roughly constant sensitivity
across time in order to treat all time steps equally in the RNTK input
comparison.

Figure 3 plots the normalized sensitivity $s(t)/{\rm max}_{t}(s(t))$ for two
data sequences of the same length $T=100$, with $s(t)$ computed numerically
for ${\bm{x}}_{t},{\bm{x}}^{\prime}_{t}\sim\mathcal{N}(0,1)$. We repeated the
experiments 10000 times; the mean of the sensitivity is shown in Figure 3.
Each of the plots shows the changes of parameters $\mathcal{S}_{\rm
ReLU}=\\{\sqrt{2},1,0,0\\}$ for $\phi=\rm ReLU$ and $\mathcal{S}_{\rm
erf}=\\{1,0.01,0.05,0\\}$ for $\phi=\rm erf$.

From Figure 3 we first observe that both $\mathrm{ReLU}$ and $\mathrm{erf}$
show similar per time step sensitivity measure $s(t)$ behavior around the
hyperparameters $\mathcal{S}_{\rm ReLU}$ and $\mathcal{S}_{\rm erf}$. If one
varies any of the weight variance parameters, the sensitivity exhibits a wide
range of behavior, and in particular with $\rm erf$. We observe that
$\sigma_{w}$ has a major influence on $s(t)$. For $\rm ReLU$, a small
decrease/increase in $\sigma_{w}$ can lead to over-sensitivity of the RNTK to
data at the last/first times steps, whereas for $\rm erf$, any changes in
$\sigma_{w}$ leads to over-sensitivity to the last time steps.

Another notable observation is the importance of $\sigma_{h}$, which is
usually set to zero for RNNs. [34] showed that a non-zero $\sigma_{h}$ acts as
a regularization that improves the performance of RNNs with the $\rm ReLU$
nonlinearity. From the sensitivity perspective, a non-zero $\sigma_{h}$
results in reducing the importance of the first time steps of the input. We
also see the same behavior in $\rm erf$, but with stronger changes as
$\sigma_{h}$ increases. Hence whenever one aims at reinforcing the input
pairwise comparisons, such parameters should be favored.

This sensitivity analysis provides a practical tool for RNTK hyperparameter
tuning. In the absence of knowledge about the data, hyperparameters should be
chosen to produce the least time varying sensitivity. If given a priori
knowledge, hyperparameters can be selected that direct the RNTK to the desired
time-steps.

##  4 Experiments

We now empirically validate the performance of the RNTK compared to classic
kernels, NTKs, and trained RNNs on both classification and regression tasks
using a large number of time series data sets. Of particular interest is the
capability of the RNTK to offer high performance even on inputs of different
lengths.

Table 1: Summary of time series classification results on 56 real-world data
sets. The RNTK outperforms classical kernels, the NTK, and trained RNNs across
all metrics. See Appendix A for detailed description of the metrics.

RNTK NTK RBF Polynomial Gaussian RNN Identity RNN GRU Acc. mean $\uparrow$
80.15% $\pm$ 15.99% 77.74% $\pm$ 16.61% 78.15% $\pm$ 16.59% 77.69% $\pm$
16.40% 55.98% $\pm$ 26.42% 63.08% $\pm$ 19.02 % 69.50% $\pm$ 22.67 P90
$\uparrow$ 92.86% 85.71% 87.60% 82.14% 28.57% 42.86% 60.71% P95 $\uparrow$
80.36% 66.07% 75.00% 67.86% 17.86% 21.43% 46.43% PMA $\uparrow$ 97.23% 94.30%
94.86% 94.23% 67.06% 78.22% 84.31% Friedman Rank $\downarrow$ 2.38 3.00 2.89
3.46 5.86 5.21 4.21

Time Series Classification. The first set of experiments considers time series
inputs of the same lengths from 56 datasets in the UCR time-series
classification data repository [13]. We restrict ourselves to selected data
sets with fewer than $1000$ training samples and fewer than $1000$ time steps
($T$) as kernel methods become rapidly intractable for larger datasets. We
compare the RNTK with a variety of other kernels, including the Radial Basis
Kernel (RBF), polynomial kernel, and NTK [24], as well as finite RNNs with
Gaussian, identity [26] initialization, and GRU [11]. We use $\phi=\rm ReLU$
for both the RNTKs and NTKs. For each kernel, we train a C-SVM [10]
classifier, and for each finite RNN we use gradient descent training. For
model hyperparameter tuning, we use 10-fold cross-validation. Details on the
data sets and experimental setup are available in Appendix 2.

We summarize the classification results over all 56 datasets in Table 1;
detailed results on each data set is available in Appendix A.2. We see that
the RNTK outperforms not only the classical kernels but also the NTK and
trained RNNs in all metrics. The results demonstrate the ability of RNTK to
provide increased performances compare to various other methods (kernels and
RNNs). The superior performance of RNTK compared to other kernels, including
NTK, can be explained by the internal recurrent mechanism present in RNTK,
allowing time-series adapted sample comparison. In addition, RNTK also
outperforms RNN and GRU. As the datasets we consider are relative small in
size, finite RNNs and GRUs that typically require large amount of data to
succeed do not perform well in our setting. An interesting future direction
would be to compare RNTK to RNN/GRU on larger datasets.

SNR

(a) Sinusoidal signal

$T_{\rm var}$

(b) Sinusoidal signal

$\sigma_{n}$

(c) Google stock value

$T_{\rm var}$

(d) Google stock value

$N_{\rm train}$

Figure 4: Performance of the RNTK on the synthetic sinusoid and real-world
Google stock price data sets compared to three other kernels. We vary the
input lengths (a,c), the input noise level (b), and training set size (d). We
compute the average SNR by repeating each experiment 1000 times. The RNTK
clearly outperforms all of the other kernels under consideration. Figure 4b
suggests that the RNTK performs better when input noise level is low
demonstrating one case where time recurrence from RNTK might be sub-optimal as
it collects and accumulate the high noise from each time step as opposed to
other kernels treating each independently.

Time Series Regression. We now validate the performance of the RNTK on time
series inputs of different lengths on both synthetic data and real data. For
both scenarios, the target is to predict the next time-step observation of the
randomly extracted windows of different length using kernel ridge regression.

We compare the RNTK to other kernels, the RBF and polynomial kernels and the
NTK. We also compare our results with a data independent predictor that
requires no training, that is simply to predict the next time step with
previous time step (PTS).

For the synthetic data experiment, we simulate 1000 samples of one period of a
sinusoid and add white Gaussian noise with default $\sigma_{n}=0.05$. From
this fixed data, we extract training set size $N_{\rm train}=20$ segments of
uniform random lengths in the range of $[T_{\rm fixed},T_{\rm fixed}+T_{\rm
var}]$ with $T_{\rm fixed}=10$. We use standard kernel ridge regression for
this task. The test set is comprised of $N_{\rm test}=5000$ obtained from
other randomly extracted segments, again of varying lengths. For the real
data, we use $975$ days of the Google stock value in the years 2014–2018. As
in the simulated signal setup above, we extract $N_{\rm train}$ segments of
different lengths from the first 700 days and test on the $N_{\rm test}$
segments from days 701 to 975. Details of the experiment are available in
Appendix A.2.

We report the predicted signal-to-noise ratio (SNR) for both datasets in
Figures 4a and 4c for various values of $T_{\rm var}$. We vary the noise level
and training set size for fixed $T_{\rm var}=10$ in Figures 4b and 4d. As we
see from Figures 4a and 4c, the RNTK offers substantial performance gains
compared to the other kernels, due to its ability to naturally deal with
variable length inputs. Moreover, the performance gap increases with the
amount of length variation of the inputs $T_{\rm var}$. Figure 4d demonstrates
that, unlike the other methods, the RNTK maintains its performance even when
the training set is small. Finally, Figure 4c demonstrates that the impact of
noise in the data on the regression performance is roughly the same for all
models but becomes more important for RNTK with a large $\sigma_{n}$; this
might be attributed to the recurrent structure of the model allowing for a
time propagation and amplification of the noise for very low SNR. These
experiments demonstrate the distinctive advantages of the RNTK over classical
kernels, and NTKs for input data sequences of varying lengths.

In the case of PTS, we expect the predictor to outperform kernel methods when
learning from the training samples is hard, due to noise in the data or small
training size which can lead to over fitting. In Figure 4a RNTK and Polynomial
kernels outperforms PTS for all values of $T_{\rm var}$, but for larger
$T_{\rm var}$, NTK and RBF under perform PTS due to the increasing detrimental
effect of zero padding. For the Google stock value, we see a superior
performance of PTS with respect to all other kernel methods due to the nature
of those data heavily relying on close past data. However, RNTK is able to
reduce the effect of over-fitting, and provide the closest results to PTS
among all kernel methods we employed, with increasing performance as the
number of training samples increase.

##  5 Conclusions

In this paper, we have derived the RNTK based on the architecture of a simple
RNN. We have proved that, at initialization, after training, and without
weight sharing, any simple RNN converges to the same RNTK. This convergence
provides new insights into the behavior of infinite-width RNNs, including how
they process different-length inputs, their training dynamics, and the
sensitivity of their output at every time step to different nonlinearities and
initializations. We have highlighted the RNTK’s practical utility by
demonstrating its superior performance on time series classification and
regression compared to a range of classical kernels, the NTK, and trained
RNNs. There are many avenues for future research, including developing RNTKs
for gated RNNs such as the LSTM [21] and investigating which of our
theoretical insights extend to finite RNNs.

## Acknowledgments

This work was supported by NSF grants CCF-1911094, IIS-1838177, and
IIS-1730574; ONR grants N00014-18-12571, N00014-20-1-2787, and
N00014-20-1-2534; AFOSR grant FA9550-18-1-0478; and a Vannevar Bush Faculty
Fellowship, ONR grant N00014-18-1-2047.

## References

  * [1] Zeyuan Allen-Zhu, Yuanzhi Li, and Yingyu Liang.  Learning and generalization in overparameterized neural networks, going beyond two layers.  In In Advances in Neural Information Processing Systems, pages 6155–6166, 2019.
  * [2] Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song.  On the convergence rate of training recurrent neural networks.  In Advances in Neural Information Processing Systems, pages 6676–6688, 2019.
  * [3] Martin Arjovsky, Amar Shah, and Yoshua Bengio.  Unitary evolution recurrent neural networks.  In International Conference on Machine Learning, pages 1120–1128, 2016.
  * [4] Sanjeev Arora, Simon S. Du, Wei Hu, Zhiyuan Li, Russ R Salakhutdinov, and Ruosong Wang.  On exact computation with an infinitely wide neural net.  In Advances in Neural Information Processing Systems, pages 8141–8150. 2019.
  * [5] Sanjeev Arora, Simon S. Du, Wei Hu, Zhiyuan Li, and Ruosong Wang.  Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks.  arXiv preprint arXiv:1901.08584, 2019.
  * [6] Sanjeev Arora, Simon S. Du, Zhiyuan Li, Ruslan Salakhutdinov, Ruosong Wang, and Dingli Yu.  Harnessing the power of infinitely wide deep nets on small-data tasks.  In International Conference on Learning Representations, 2020.
  * [7] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  Proceedings of the National Academy of Sciences, 116(32):15849–15854, 2019.
  * [8] Yoshua. Bengio, Patrice. Simard, and Paolo. Frasconi.  Learning long-term dependencies with gradient descent is difficult.  IEEE Trans. Neural Networks, 5(2):157–166, 1994.
  * [9] Erwin Bolthausen.  An iterative construction of solutions of the tap equations for the sherrington–kirkpatrick model.  Communications in Mathematical Physics, 325(1):333–366, 2014.
  * [10] Chih-Chung Chang and Chih-Jen Lin.  LIBSVM: A library for support vector machines.  ACM Transactions on Intelligent Systems and Technology, 2:27:1–27:27, 2011.  Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm.
  * [11] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  arXiv preprint arXiv:1406.1078, 2014.
  * [12] Youngmin Cho and Lawrence K Saul.  Kernel methods for deep learning.  In Advances in Neural Information Processing Systems, pages 342–350, 2009.
  * [13] Hoang Anh Dau, Eamonn Keogh, Kaveh Kamgar, Chin-Chia Michael Yeh, Yan Zhu, Shaghayegh Gharghabi, Chotirat Ann Ratanamahatana, Yanping Chen, Bing Hu, Nurjahan Begum, Anthony Bagnall, Abdullah Mueen, Gustavo Batista, and ML Hexagon.  The UCR time series classification archive, 2019.
  * [14] Simon S. Du, Kangcheng Hou, Russ R Salakhutdinov, Barnabas Poczos, Ruosong Wang, and Keyulu Xu.  Graph neural tangent kernel: Fusing graph neural networks with graph kernels.  In Advances in Neural Information Processing Systems, pages 5724–5734, 2019.
  * [15] Simon S. Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai.  Gradient descent finds global minima of deep neural networks.  In International Conference on Machine Learning, pages 1675–1685, 2019.
  * [16] David Duvenaud, Oren Rippel, Ryan Adams, and Zoubin Ghahramani.  Avoiding pathologies in very deep networks.  In Artificial Intelligence and Statistics, pages 202–210, 2014\.
  * [17] Jeffrey L. Elman.  Finding structure in time.  Cognitive Science, 14(2):179–211, 1990.
  * [18] M. Fernández-Delgado, M.S. Sirsat, E. Cernadas, S. Alawadi, S. Barro, and M. Febrero-Bande.  An extensive experimental survey of regression methods.  Neural Networks, 111:11 – 34, 2019.
  * [19] Adrià Garriga-Alonso, Carl Edward Rasmussen, and Laurence Aitchison.  Deep convolutional networks as shallow gaussian processes.  In International Conference on Learning Representations, 2019.
  * [20] Mikael Henaff, Arthur Szlam, and Yann LeCun.  Recurrent orthogonal networks and long-memory tasks.  In International Conference on Machine Learning, pages 2034–2042. PMLR, 2016.
  * [21] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural Computation, 9(8):1735–1780, 1997.
  * [22] Wei Hu, Zhiyuan Li, and Dingli Yu.  Simple and effective regularization methods for training on noisily labeled data with generalization guarantee.  In International Conference on Learning Representations, 2020.
  * [23] Kaixuan Huang, Yuqing Wang, Molei Tao, and Tuo Zhao.  Why do deep residual networks generalize better than deep feedforward networks?–a neural tangent kernel perspective.  arXiv preprint arXiv:2002.06262, 2020.
  * [24] Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  In Advances in neural information processing systems, pages 8571–8580, 2018.
  * [25] Li Jing, Yichen Shen, Tena Dubcek, John Peurifoy, Scott Skirlo, Yann LeCun, Max Tegmark, and Marin Soljačić.  Tunable efficient unitary neural networks (eunn) and their application to rnns.  In International Conference on Machine Learning, pages 1733–1741, 2017.
  * [26] Quoc V Le, Navdeep Jaitly, and Geoffrey E Hinton.  A simple way to initialize recurrent networks of rectified linear units.  arXiv preprint arXiv:1504.00941, 2015.
  * [27] Jaehoon Lee, Samuel S Schoenholz, Jeffrey Pennington, Ben Adlam, Lechao Xiao, Roman Novak, and Jascha Sohl-Dickstein.  Finite versus infinite neural networks: an empirical study.  arXiv preprint arXiv:2007.15801, 2020.
  * [28] Jaehoon Lee, Jascha Sohl-dickstein, Jeffrey Pennington, Roman Novak, Sam Schoenholz, and Yasaman Bahri.  Deep neural networks as gaussian processes.  In International Conference on Learning Representations, 2018.
  * [29] Jaehoon Lee, Lechao Xiao, Samuel Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, and Jeffrey Pennington.  Wide neural networks of any depth evolve as linear models under gradient descent.  In Advances in Neural Information Processing Systems, pages 8572–8583, 2019.
  * [30] Radford M Neal.  Bayesian Learning for Neural Networks. PhD thesis, University of Toronto, 1995.
  * [31] Behnam Neyshabur, Zhiyuan Li, Srinadh Bhojanapalli, Yann LeCun, and Nathan Srebro.  The role of over-parametrization in generalization of neural networks.  In International Conference on Learning Representations, 2019.
  * [32] Roman Novak, Yasaman Bahri, Daniel A. Abolafia, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Sensitivity and generalization in neural networks: an empirical study.  In International Conference on Learning Representations, 2018.
  * [33] Roman Novak, Lechao Xiao, Yasaman Bahri, Jaehoon Lee, Greg Yang, Daniel A. Abolafia, Jeffrey Pennington, and Jascha Sohl-dickstein.  Bayesian deep convolutional networks with many channels are gaussian processes.  In International Conference on Learning Representations, 2019.
  * [34] Zichao Wang, Randall Balestriero, and Richard Baraniuk.  A max-affine spline perspective of recurrent neural networks.  In International Conference on Learning Representations, 2018.
  * [35] Greg Yang.  Scaling limits of wide neural networks with weight sharing: Gaussian process behavior, gradient independence, and neural tangent kernel derivation.  arXiv preprint arXiv:1902.04760, 2019.
  * [36] Greg Yang.  Tensor programs I: Wide feedforward or recurrent neural networks of any architecture are gaussian processes.  arXiv preprint arXiv:1910.12478, 2019.
  * [37] Greg Yang.  Tensor programs II: Neural tangent kernel for any architecture.  arXiv preprint arXiv:2006.14548, 2020.
  * [38] Greg Yang.  Tensor programs III: Neural matrix laws.  arXiv preprint arXiv:2009.10685, 2020.
  * [39] Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu.  Stochastic gradient descent optimizes over-parameterized deep ReLU networks.  arXiv preprint arXiv:1811.08888, 2018.

##  Appendix A Experiment Details

###  A.1 Time series classification

Table 2: Test accuracy of each model on 56 time series data set from UCR time-series classification data repository [13]. Dataset | RNTK | NTK | RBF | POLY | Gaussian RNN | Identity RNN | GRU
---|---|---|---|---|---|---|---
Strawberry | 98.38 | 97.57 | 97.03 | 96.76 | 94.32 | 75.4 | 91.62
ProximalPhalanxOutlineCorrect | 89 | 87.97 | 87.29 | 86.94 | 82.81 | 74.57 | 86.94
PowerCons | 97.22 | 97.22 | 96.67 | 91.67 | 96.11 | 95 | 99.44
Ham | 70.48 | 71.63 | 66.67 | 71.43 | 53.33 | 60 | 60.95
SmallKitchenAppliances | 67.47 | 38.4 | 40.27 | 37.87 | 60.22 | 76 | 71.46
ScreenType | 41.6 | 43.2 | 43.47 | 38.4 | 40 | 41.06 | 36.26
MiddlePhalanxOutlineCorrect | 57.14 | 57.14 | 48.7 | 64.29 | 76.28 | 57.04 | 74.57
RefrigerationDevices | 46.93 | 37.07 | 36.53 | 41.07 | 36 | 50.93 | 46.66
Yoga | 84.93 | 84.63 | 84.63 | 84.87 | 46.43 | 76.66 | 61.83
Computers | 59.2 | 55.2 | 58.8 | 56.4 | 53.2 | 55.2 | 58.8
ECG5000 | 93.76 | 94.04 | 93.69 | 93.96 | 88.4 | 93.15 | 93.26
Fish | 90.29 | 84 | 85.71 | 88 | 28 | 38.28 | 24
UWaveGestureLibraryX | 79.59 | 78.7 | 78.48 | 65.83 | 55.97 | 75.34 | 73.64
UWaveGestureLibraryY | 71.56 | 70.63 | 70.35 | 70.32 | 44.5 | 65.18 | 65.38
UWaveGestureLibraryZ | 73.95 | 73.87 | 72.89 | 71.94 | 43.29 | 67.81 | 70.32
StarLightCurves | 95.94 | 96.19 | 94.62 | 94.44 | 82.13 | 86.81 | 96.15
CricketX | 60.51 | 59.49 | 62.05 | 62.56 | 8.46 | 63.58 | 26.41
CricketY | 63.85 | 58.97 | 60.51 | 59.74 | 15.89 | 59.23 | 36.15
CricketZ | 60.26 | 59.23 | 62.05 | 59.23 | 8.46 | 57.94 | 41.28
DistalPhalanxOutlineCorrect | 77.54 | 77.54 | 75.36 | 73.91 | 69.92 | 69.56 | 75
Worms | 57.14 | 50.65 | 55.84 | 50.65 | 35.06 | 49.35 | 41.55
SyntheticControl | 98.67 | 96.67 | 98 | 97.67 | 92.66 | 97.66 | 99
Herring | 56.65 | 59.38 | 59.38 | 59.38 | 23.28 | 59.37 | 59.38
MedicalImages | 74.47 | 73.29 | 75.26 | 74.61 | 48.15 | 64.86 | 69.07
SwedishLeaf | 90.56 | 91.04 | 91.36 | 90.72 | 59.2 | 45.92 | 91.04
ChlorineConcentration | 90.76 | 77.27 | 86.35 | 91.54 | 65.99 | 55.75 | 61.14
SmoothSubspace | 96 | 87.33 | 92 | 86.67 | 94 | 95.33 | 92.66
TwoPatterns | 94.25 | 90.45 | 91.25 | 93.88 | 99.7 | 99.9 | 100
Faceall | 74.14 | 83.33 | 83.25 | 82.43 | 53.66 | 70.53 | 70.65
DistalPhalanxTW | 66.19 | 69.78 | 66.91 | 67.37 | 67.62 | 64.74 | 69.06
MiddlePhalanxTW | 57.79 | 61.04 | 59.74 | 60.39 | 58.44 | 58.44 | 59.09
FacesUCR | 81.66 | 80.2 | 80.34 | 82.98 | 53.21 | 75.26 | 79.46
OliveOil | 90 | 86.67 | 86.67 | 83.33 | 66.66 | 40 | 40
UMD | 91.67 | 92.36 | 97.22 | 90.97 | 44.44 | 71.52 | 100
nsectEPGRegular | 99.6 | 99.2 | 99.6 | 96.79 | 100 | 100 | 98.39
Meat | 93.33 | 93.33 | 93.33 | 93.33 | 0.55 | 55 | 33.33
Lightning2 | 78.69 | 73.77 | 70.49 | 68.85 | 45.9 | 70.49 | 67.21
Lightning7 | 61.64 | 60.27 | 63.01 | 60.27 | 23.28 | 69.86 | 76.71
Car | 83.33 | 78.83 | 80 | 80 | 23.33 | 58.33 | 26.66
GunPoint | 98 | 95.33 | 95.33 | 94 | 82 | 74.66 | 80.66
Arrowhead | 80.57 | 83.43 | 80.57 | 74.86 | 48 | 56 | 37.71
Coffee | 100 | 100 | 92.86 | 92.86 | 100 | 42.85 | 57.14
Trace | 96 | 81 | 76 | 76 | 70 | 71 | 100
ECG200 | 93 | 89 | 89 | 86 | 86 | 72 | 76
plane | 98.1 | 96.19 | 97.14 | 97.14 | 96.19 | 84.76 | 96.19
GunPointOldVersusYoung | 98.73 | 97.46 | 98.73 | 94.6 | 53.96 | 52.38 | 98.41
GunPointMaleVersusFemale | 99.05 | 99.68 | 99.37 | 99.68 | 68.67 | 52.53 | 99.68
GunPointAgeSpan | 96.52 | 94.62 | 95.89 | 93.99 | 47.78 | 47.78 | 95.56
FreezerRegularTrain | 97.44 | 94.35 | 96.46 | 96.84 | 76.07 | 7.5 | 86.59
SemgHandSubjectCh2 | 84.22 | 85.33 | 86.14 | 86.67 | 20 | 36.66 | 89.11
WormsTwoClass | 62.34 | 62.34 | 61.04 | 59.74 | 51.94 | 46.75 | 57.14
Earthquakes | 74.82 | 74.82 | 74.82 | 74.82 | 65.46 | 76.97 | 76.97
FiftyWords | 68.57 | 68.57 | 69.67 | 68.79 | 34.28 | 60.21 | 65.27
Beef | 90 | 73.33 | 83.33 | 93.33 | 26.67 | 46.67 | 36.67
Adiac | 76.63 | 71.87 | 73.40 | 77.75 | 51.4 | 16.88 | 60.61
WordSynonyms | 57.99 | 58.46 | 61.13 | 62.07 | 17.71 | 45.77 | 53.76

#### Kernel methods settings.

We used RNTK, RBF, polynomial and NTK [24]. For data pre-processing, we
normalized the norm of each ${\bm{x}}$ to 1. For training we used C-SVM in
LIBSVM library [10] and for hyperparameter selection we performed 10-fold
validation for splitting the training data into 90% training set and 10%
validation test. We then choose the best performing set of hyperparameters on
all the validation sets, retrain the models with the best set of
hyperparameters on the entire training data and finally report the performance
on the unseen test data. The performance of all kernels on each data set is
shown in table 2.

For C-SVM we chose the cost function value

$\displaystyle C\in\\{0.01,0.1,1,10,100\\}$

and for each kernel we used the following hyperparameter sets

  * •

RNTK: We only used single layer RNTK, we $\phi={\rm ReLU}$ and the following
hyperparameter sets for the variances:

$\displaystyle\sigma_{w}\in\\{1.34,1.35,1.36,1.37,1.38,1.39,1.40,1.41,1.42,\sqrt{2},1.43,1.44,1.45,1.46,1.47\\}$
$\displaystyle\sigma_{u}=1$
$\displaystyle\sigma_{b}\in\\{0,0.01,0.05,0.1,0.2,0.3,0.4,0.5,0.7,0.9,1,2\\}$
$\displaystyle\sigma_{h}\in\\{0,0.01,0.1,0.5,1\\}$

  * •

NTK: The formula for NTK of $L$-layer MLP [24] for
${\bm{x}},{\bm{x}}^{\prime}\in\mathbb{R}^{m}$ is:

$\displaystyle\Sigma^{(1)}$
$\displaystyle=\frac{\sigma_{w}^{2}}{m}\langle{\bm{x}},{\bm{x}}^{\prime}\rangle+\sigma_{b}^{2}$
$\displaystyle\Sigma^{(\ell)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi}[K^{(\ell)}({\bm{x}},{\bm{x}}^{\prime})]+\sigma_{b}^{2}$
$\displaystyle\ell\in[L]$
$\displaystyle\dot{\Sigma}^{(\ell)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi^{\prime}}[K^{(\ell+1)}({\bm{x}},{\bm{x}}^{\prime})]$
$\displaystyle\ell\in[L]$
$\displaystyle{\bm{K}}^{(\ell)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\left[{\begin{array}[]{cc}\Sigma^{(\ell-1)}({\bm{x}},{\bm{x}})&\Sigma^{(\ell-1)}({\bm{x}},{\bm{x}}^{\prime})\\\
\Sigma^{(\ell-1)}({\bm{x}},{\bm{x}}^{\prime})&\Sigma^{(\ell-1)}({\bm{x}}^{\prime},{\bm{x}}^{\prime})\\\
\end{array}}\right]$ (25)
$\displaystyle\mathcal{K}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi}[K^{(L+1)}({\bm{x}},{\bm{x}}^{\prime})]$
$\displaystyle k_{\rm NTK}$
$\displaystyle=\sum_{\ell=1}^{L}\left(\Sigma^{(\ell)}({\bm{x}},{\bm{x}}^{\prime})\prod_{\ell^{\prime}=\ell}^{L}\dot{\Sigma}^{(\ell)}({\bm{x}},{\bm{x}}^{\prime})\right)+\mathcal{K}({\bm{x}},{\bm{x}}^{\prime})$

and we used the following hyperparamters

$\displaystyle L\in[10]$
$\displaystyle\sigma_{w}\in\\{0.5,1,\sqrt{2},2,2.5,3\\}$
$\displaystyle\sigma_{b}\in\\{0,0.01,0.1,0.2,0.5,0.8,1,2,5\\}$

  * •

RBF:

$\displaystyle k_{\rm
RBF}({\bm{x}},{\bm{x}}^{\prime})=e^{(-\alpha\|{\bm{x}}-{\bm{x}}^{\prime}\|^{2}_{2})}$
$\displaystyle\alpha\in\\{0.01,0.05,0.1,0.2,0.5,0.6,0.7,0.8,1,2,3,4,5,10,20,30,40,100\\}$

  * •

Polynomial:

$\displaystyle k_{\rm
Polynomial}({\bm{x}},{\bm{x}}^{\prime})=(r+\langle{\bm{x}},{\bm{x}}^{\prime}\rangle)^{d}$
$\displaystyle d\in[5]$ $\displaystyle r\in\\{0,0.1,0.2,0.5,1,2\\}$

#### Finite-width RNN settings.

We used 3 different RNNs. The first is a ReLU RNN with Gaussian initialization
with the same NTK initialization scheme, where parameter variances are
$\sigma_{w}=\sigma_{v}=\sqrt{2}$, $\sigma_{u}=1$ and $\sigma_{b}=0$. The
second is a ReLU RNN with identity initialization following [26]. The third is
a GRU [11] with uniform initialization. All models are trained with RMSProp
algorithm for 200 epochs. Early stopping is implemented when the validation
set accuracy does not improve for 5 consecutive epochs.

We perform standard 5-fold cross validation. For each RNN architecture we used
hyperparamters of number of layer, number of hidden units and learning rate as

$\displaystyle L\in\\{1,2\\}$ $\displaystyle n\in\\{50,100,200,500\\}$
$\displaystyle\eta\in\\{0.01,0.001,0.0001,0.00001\\}$

#### Metrics descriptions

First, only in this paragraph, let $i\in\\{1,2,...,N\\}$ index a total of $N$
datasets and $j\in\\{1,2,...,M\\}$ index a total of $M$ classifiers. Let
$y_{ij}$ be the accuracy of the $j$-th classifer on the $i$-th dataset. We
reported results on 4 metrics: average accuracy (Acc. mean), P90, P95, PMA and
Friedman Rank. P90 and P95 is the fraction of datasets that the classifier
achieves at least 90% and 95% of the maximum achievable accuracy for each
dataset, i.e.,

$\displaystyle P90_{j}=\frac{1}{N}\sum_{i}\mathbf{1}(y_{ij}\geq
0.9(\underset{j}{\rm max}\,y_{ij}))\,.$ (26)

PMA is the accuracy of the classifier on a dataset divided by the maximum
achievable accuracy on that dataset, averaged over all datasets:

$\displaystyle{\rm PMA}_{j}=\frac{1}{N}\sum_{i}\frac{y_{ij}}{\underset{j}{\rm
max}\,y_{ij}}\,.$ (27)

Friedman Rank [18] first ranks the accuracy of each classifier on each dataset
and then takes the average of the ranks for each classifier over all datasets,
i.e.,

$\displaystyle{\rm FR}_{j}=\frac{1}{N}\sum_{i}r_{ij}\,,$ (28)

where $r_{ij}$ is the ranking of the $j$-th classifier on the $i$-th dataset.

Note that a better classifier achieves a lower Friedman Rank, Higher P/90 and
PMA.

Remark. In order to provide insight into the performance of RNTK in long time
steps setting, we picked two datasets with more that 1000 times steps:
SemgHandSubjectCh2 ($T=1024$) and StarLightCurves ($T=1024$).

###  A.2 Time Series Regression

For time series regression, we used the 5-fold validation of training set and
same hyperparamter sets for all kernels. For training we kernel ridge
regression with ridge term chosen form

$\displaystyle\lambda\in\\{0,0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.8,1,2,3,4,5,6,7,8,10,100\\}$

##  Appendix B Proofs for Theorems 1 and 3: RNTK Convergence at Initialization

###  B.1 Preliminary: Netsor programs

Calculation of NTK in any architecture relies on finding the GP kernels that
correspond to each pre-activation and gradient layers at initialization. For
feedforward neural networks with $n_{1},\dots,n_{L}$ number of neurons
(channels in CNNs) at each layer the form of this GP kernels can be calculated
via taking the limit of $n_{1},\dots,n_{L}$ sequentially one by one. The proof
is given by induction, where by conditioning on the previous layers, each
entry of the current layer is sum of infinite i.i.d Gaussian random variables,
and based on Central Limit Theorem (CLT), it becomes a Gaussian process with
kernel calculated based on the previous layers. Since the first layer is an
affine transformation of input with Gaussian weights, it is a Gaussian process
and the proof is completed. See [28, 16, 33, 19] for a formal treatment.
However, due to weight-sharing, sequential limit is not possible and condoning
on previous layers does not result in i.i.d. weights. Hence the aforementioned
arguments break. To deal with it, in [35] a proof using Gaussian conditioning
trick [9] is presented which allows use of recurrent weights in a network.
More precisely, it has been demonstrated than neural networks (without batch
normalization) can be expressed and a series of matrix multiplication and
(piece wise) nonlinearity application, generally referred as _Netsor
programs_. It has been shown that any architecture that can be expressed as
_Netsor programs_ that converge to GPs as width goes to infinity in the same
rate, which a general rule to obtain the GP kernels. For completeness of this
paper, we briefly restate the results from [35] which we will use later for
calculation derivation of RNTK.

There are 3 types of variables in _Netsor programs_ ; $A$-vars, $G$-vars and
$H$-vars. $A$-vars are matrices and vectors with i.i.d Gaussian entries,
$G$-vars are vectors introduced by multiplication of a vector by an $A$-var
and $H$-vars are vectors after coordinate wise nonlinearities is applied to
$G$-vars. Generally, $G$-vars can be thought of as pre-activation layers which
are asymptotically treated as a Gaussian distributed vectors, $H$-vars as
after-activation layers and $A$-vars are the weights. Since in neural networks
inputs are immediately multiplied by a weight matrix, it can be thought of as
an $G$-var, namely ${\bm{g}}_{in}$. Generally _Netsor programs_ supports
$G$-vars with different dimension, however the asymptotic behavior of a neural
networks described by _Netsor programs_ does not change under this degree of
freedom, as long as they go to infinity at the same rate. For simplicity, let
the $G$-vars and $H$-vars have the same dimension $n$ since the network of
interest is RNN and all pre-activation layers have the same dimension. We
introduce the _Netsor programs_ under this simplification. To produce the
output of a neural network, _Netsor programs_ receive a set of $G$-vars and
$A$-vars as input, and new variables are produced sequentially using the three
following operators:

  * •

Matmul : multiplication of an $A$-var: ${\bm{A}}$ with an $H$-var: ${\bm{h}}$,
which produce a new $G$-var, ${\bm{g}}$.

$\displaystyle{\bm{g}}={\bm{A}}{\bm{h}}$ (29)

  * •

Lincomp: Linear combination of $G$-vars, ${\bm{g}}^{i},\,\,1\leq i\leq k$ ,
with coefficients $a^{i}\in\mathbb{R}\,\,1\leq i\leq k$ which produce of new
$G$-var:

$\displaystyle{\bm{g}}=\sum_{i=1}^{k}a^{i}{\bm{g}}^{i}$ (30)

  * •

Nonlin: creating a new $H$-var, ${\mathbf{h}}$, by using a nonlinear function
$\phi:\mathbb{R}^{k}\rightarrow\mathbb{R}$ that act coordinate wise on a set
of $G$-vars, ${\bm{g}}^{i},\,\,1\leq i\leq k$ :

$\displaystyle{\bm{h}}=\varphi({\bm{g}}^{1},\dots,{\bm{g}}^{k})$ (31)

Any output of the neural network $y\in\mathbb{R}$ should be expressed as inner
product of a new $A$-var which has not been used anywhere else in previous
computations and an $H$-var:

$\displaystyle y={\bm{v}}^{\top}{\bm{h}}$

Any other output can be produced by another ${\bm{v}}^{\prime}$ and
${\bm{h}}^{\prime}$ (possibility the same ${\bm{h}}$ or ${\bm{v}}$).

It is assumed that each entry of any $A$-var : ${\bm{A}}\in\mathbb{R}^{n\times
n}$ in the _netsor programs_ computations is drawn from
$\mathcal{N}(0,\frac{\sigma_{a}^{2}}{n})$ and the input $G$-vars are Gaussian
distributed. The collection of a specific entry of all $G$-vars of in the
_netsor program_ converges in probability to a Gaussian vector
$\\{[{\bm{g}}^{1}]_{i},\dots,[{\bm{g}}^{k}]_{i}\\}\sim\mathcal{N}({\bm{\mu}},{\bm{\Sigma}})$
for all $i\in[n]$ as $n$ goes to infinity.
Let $\mu({\bm{g}}):={\mathbb{E}}\big{[}[{\bm{g}}]_{i}\big{]}$ be the mean of a
$G$-var and
$\Sigma({\bm{g}},{\bm{g}}^{\prime}):={\mathbb{E}}\big{[}[{\bm{g}}]_{i}\bm{\cdot}[{\bm{g}}^{\prime}]_{i}\big{]}$
be the covariance between any two $G$-vars. The general rule for
$\mu({\bm{g}})$ is given by the following equations:

$\displaystyle\mu({\bm{g}})=\begin{dcases}\mu^{in}({\bm{g}})&\text{if
}{\bm{g}}\text{ is input}\\\ \sum_{i=1}^{k}a^{i}\mu({\bm{g}}^{i})&\text{if
}{\bm{g}}=\sum_{i=1}^{k}a^{i}{\bm{g}}^{i}\\\ 0&\text{otherwise}\\\
\end{dcases}$ (32)

For ${\bm{g}}$ and ${\bm{g}}^{\prime}$, let
$\mathcal{G}=\\{{\bm{g}}^{1},\dots,{\bm{g}}^{r}\\}$ be the set of $G$-vars
that has been introduced _before_ ${\bm{g}}$ and ${\bm{g}}^{\prime}$ with
distribution
$\mathcal{N}({\bm{\mu}}_{\mathcal{G}},{\bm{\Sigma}}_{\mathcal{G}})$, where
${\bm{\Sigma}}_{\mathcal{G}}\in\mathbb{R}^{|\mathcal{G}|\times|\mathcal{G}|}$
containing the pairwise covariances between the $G$-vars.
$\Sigma({\bm{g}},{\bm{g}}^{\prime})$ is calculated via the following rules:

$\displaystyle\Sigma({\bm{g}},{\bm{g}}^{\prime})=\begin{dcases}\Sigma^{\rm
in}({\bm{g}},{\bm{g}}^{\prime})&\text{if }{\bm{g}}\text{ and
}{\bm{g}}^{\prime}\text{ are inputs}\\\
\sum_{i=1}^{k}a^{i}\Sigma({\bm{g}}^{i},{\bm{g}}^{\prime})&\text{if
}{\bm{g}}=\sum_{i=1}^{k}a^{i}{\bm{g}}^{i}\\\
\sum_{i=1}^{k}a^{i}\Sigma({\bm{g}},{\bm{g}}^{i})&\text{if
}{\bm{g}}^{\prime}=\sum_{i=1}^{k}a^{i}{\bm{g}}^{i}\\\
\sigma^{2}_{A}\mathop{\mathbb{E}}_{{\mathbf{z}}\sim\mathcal{N}({\bm{\mu}},{\bm{\Sigma}}_{\mathcal{G}})}[\varphi({\mathbf{z}})\bar{\varphi}({\mathbf{z}})]&\text{if
}{\bm{g}}={\bm{A}}{\bm{h}}\text{ and
}{\bm{g}}^{\prime}={\bm{A}}{\bm{h}}^{\prime}\\\
0&\text{otherwise}\end{dcases}$ (33)

Where ${\bm{h}}=\varphi({\bm{g}}^{1},\dots,{\bm{g}}^{r})$ and
${\bm{h}}^{\prime}=\bar{\varphi}({\bm{g}}^{1},\dots,{\bm{g}}^{r})$ are
functions of $G$-vars in $\mathcal{G}$ from possibly different nonlinearities.
This set of rules presents a recursive method for calculating the GP kernels
in a network where the recursive formula starts from data dependent quantities
$\Sigma^{\rm in}$ and $\mu^{in}$ which are given.

All the above results holds when the nonlinearities are bounded uniformly by
$e^{(cx^{2-\alpha})}$ for some $\alpha>0$ and when their derivatives exist.

Standard vs. NTK initialization. The common practice (which _netsor programs_
uses) is to initialize DNNs weights $[{\bm{A}}]_{i,j}$ with
$\mathcal{N}(0,\frac{\sigma_{a}}{\sqrt{n}})$ (known as _standard_
initialization) where generally $n$ is the number of units in the previous
layer. In this paper we have used a different parameterization scheme as used
in [24] and we factor the standard deviation as shown in 3 and initialize
weights with standard standard Gaussian. This approach does not change the the
forward computation of DNN, but normalizes the backward computation (when
computing the gradients) by factor $\frac{1}{n}$, otherwise RNTK will be
scales by $n$. However this problem can be solved by scaling the step size by
$\frac{1}{n}$ and there is no difference between _NTK_ and _standard_
initialization [29].

###  B.2 Proof for Theorem 1: Single layer Case

We first derive the RNTK in a simpler setting, i.e., a single layer and single
output RNN. We then generalize the results to multi-layer and multi-output
RNNs. We drop the layer index $\ell$ to simplify notation. From 3 and 4, the
forward pass for computing the output under NTK initialization for each input
${\bm{x}}=\\{{\bm{x}}_{t}\\}_{t=1}^{T}$ is given by:

$\displaystyle{\bm{g}}^{\left(t\right)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{w}}{\sqrt{m}}{\mathbf{W}}{\bm{h}}^{\left(t-1\right)}({\bm{x}})+\frac{\sigma_{u}}{\sqrt{n}}{\mathbf{U}}{\bm{x}}_{t}+\sigma_{b}{\mathbf{b}}$
(34) $\displaystyle{\bm{h}}^{\left(t\right)}({\bm{x}})$
$\displaystyle=\phi\left({\bm{g}}^{\left(t\right)}({\bm{x}})\right)$ (35)
$\displaystyle f_{\theta}({\bm{x}})$
$\displaystyle=\frac{\sigma_{v}}{\sqrt{n}}{\bm{v}}^{\top}{\bm{h}}^{(T)}({\bm{x}})$
(36)

Note that (34), (35) and (36) use all the introduced operators introduced in
29, 30 and 31 given input variables
${\mathbf{W}},\\{{\mathbf{U}}{\bm{x}}_{t}\\}_{t=1}^{T},{\mathbf{b}},{\bm{v}}$
and ${\bm{h}}^{(0)}({\bm{x}})$.

First, we compute the kernels of forward pass
$\Sigma^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$ and backward pass
$\Pi^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$ introduced in (6) and (7)
for two input ${\bm{x}}$ and ${\bm{x}}^{\prime}$. Note that based on (32) the
mean of all variables is zero since the inputs are all zero mean. In the
forward pass for the intermediate layers we have:

$\displaystyle\Sigma^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})=\Sigma({\bm{g}}^{\left(t\right)}({\bm{x}}),{\bm{g}}^{(t^{\prime})}({\bm{x}}^{\prime}))$
$\displaystyle=\Sigma\left(\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{\left(t-1\right)}({\bm{x}})+\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}_{t}+\sigma_{b}{\mathbf{b}},\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t^{\prime}-1)}({\bm{x}}^{\prime})+\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}^{\prime}_{t^{\prime}}+\sigma_{b}{\mathbf{b}}\right)$
$\displaystyle=\Sigma\left(\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t-1)}({\bm{x}}),\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t^{\prime}-1)}({\bm{x}}^{\prime})\right)+\Sigma^{\rm
in}\left(\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}_{t},\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}^{\prime}_{t^{\prime}}\right)+\Sigma^{\rm
in}\left(\sigma_{b}{\bm{b}},\sigma_{b}{\bm{b}}\right).$

We have used the second and third rule in (33) to expand the formula, We have
also used the first and fifth rule to set the cross term to zero, i.e.,

$\displaystyle\Sigma\left(\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t-1)}({\bm{x}}),\frac{\sigma_{u}}{\sqrt{n}}{\mathbf{U}}{\bm{x}}^{\prime}_{t^{\prime}}\right)$
$\displaystyle=0$
$\displaystyle\Sigma\left(\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t-1)}({\bm{x}}),\sigma_{b}{\mathbf{b}}\right)$
$\displaystyle=0$
$\displaystyle\Sigma\left(\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}_{t},\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t^{\prime}-1)}({\bm{x}}^{\prime})\right)$
$\displaystyle=0$
$\displaystyle\Sigma\left(\sigma_{b}{\mathbf{b}},\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t^{\prime}-1)}({\bm{x}}^{\prime})\right)$
$\displaystyle=0$ $\displaystyle\Sigma^{\rm
in}\left(\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}_{t},\sigma_{b}{\mathbf{b}}\right)$
$\displaystyle=0$ $\displaystyle\Sigma^{\rm
in}\left(\sigma_{b}{\mathbf{b}},\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}^{\prime}_{t^{\prime}}\right)$
$\displaystyle=0.$

For the non-zero terms we have

$\displaystyle\Sigma^{\rm
in}\left(\sigma_{b}{\bm{b}},\sigma_{b}{\bm{b}}\right)$
$\displaystyle=\sigma_{b}^{2}$ $\displaystyle\Sigma^{\rm
in}\left(\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}_{t},\frac{\sigma_{u}}{\sqrt{m}}{\mathbf{U}}{\bm{x}}^{\prime}_{t^{\prime}}\right)$
$\displaystyle=\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{t^{\prime}}\rangle,$

which can be achieved by straight forward computation. If $t\neq 1$ and
$t^{\prime}\neq 1$, by using the forth rule in (33) we have

$\displaystyle\Sigma\left(\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t-1)}({\bm{x}}),\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}{\bm{h}}^{(t^{\prime}-1)}({\bm{x}}^{\prime})\right)=\sigma_{w}^{2}\mathop{\mathbb{E}}_{{\mathbf{z}}\sim\mathcal{N}(0,{\bm{K}}^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime}))}[\phi({\textnormal{z}}_{1})\phi({\textnormal{z}}_{2})]=\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\big{]}.$

With ${\bm{K}}^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$ defined in (18).
Otherwise, it will be zero by the fifth rule (if $t$ or $t$ = 1) .

Here the set of previously introduced $G$-vars is
$\mathcal{G}=\big{\\{}\\{{\bm{g}}^{(\alpha)}({\bm{x}})\\},{\mathbf{U}}{\bm{x}}_{\alpha}\\}_{\alpha=1}^{t-1},\\{{\bm{g}}^{(\alpha^{\prime})}({\bm{x}}^{\prime}),{\mathbf{U}}{\bm{x}}^{\prime}_{\alpha^{\prime}}\\}_{\alpha^{\prime}=1}^{t^{\prime}-1},{\bm{h}}^{(0)}({\bm{x}}),{\bm{h}}^{(0)}({\bm{x}}^{\prime})\big{\\}}$,
but the dependency is only on the last layer $G$-vars,
$\varphi(\\{{\bm{g}}:{\bm{g}}\in\mathcal{G}\\})=\phi({\bm{g}}^{(t-1)}({\bm{x}}))$,
$\bar{\varphi}((\\{{\bm{g}}:{\bm{g}}\in\mathcal{G}\\}))=\phi({\bm{g}}^{(t^{\prime}-1)}({\bm{x}}^{\prime}))$,
leading the calculation to the operator defined in (10). As a result

$\displaystyle\Sigma^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\big{]}+\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{t^{\prime}}\rangle+\sigma_{b}^{2}.$

To complete the recursive formula, using the same procedure for the first
layers we have

$\displaystyle\Sigma^{(1,1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\sigma_{h}^{2}1_{({\bm{x}}={\bm{x}}^{\prime})}+\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{1},{\bm{x}}^{\prime}_{1}\rangle+\sigma_{b}^{2},$
$\displaystyle\Sigma^{(1,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{1},{\bm{x}}^{\prime}_{t^{\prime}}\rangle+\sigma_{b}^{2},$
$\displaystyle\Sigma^{(t,1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{1}\rangle+\sigma_{b}^{2}.$

The output GP kernel is calculated via

$\displaystyle\mathcal{K}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(T+1,T^{\prime}+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}$

The calculation of the gradient vectors
${\bm{\delta}}^{(t)}({\bm{x}})=\sqrt{n}\big{(}\nabla_{{\bm{g}}^{(t)}({\bm{x}})}f_{\theta}({\bm{x}})\big{)}$
in the backward pass is given by

$\displaystyle{\bm{\delta}}^{(T)}({\bm{x}})$
$\displaystyle=\sigma_{v}{\bm{v}}\odot\phi^{\prime}({\bm{g}}^{(T)}({\bm{x}}))$
$\displaystyle{\bm{\delta}}^{(t)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}^{\top}\Big{(}\phi^{\prime}({\bm{g}}^{(t)}({\bm{x}}))\odot{\bm{\delta}}^{(t+1)}({\bm{x}})\Big{)}\hskip
20.0ptt\in[T-1]$

To calculate the backward pass kernels, we rely on the following Corollary
from [38]

###### Corollary 1

In infinitely wide neural networks weights used in calculation of back
propagation gradients (${\mathbf{W}}^{\top}$) is an i.i.d copy of weights used
in forward propagation (${\mathbf{W}}$) as long as the last layer weight
(${\bm{v}}$) is sampled independently from other parameters and has mean 0.

The immediate result of Corollary 1 is that ${\bm{g}}^{(t)}({\bm{x}})$ and
${\bm{\delta}}^{(t)}({\bm{x}})$ are two independent Gaussian vector as their
covariance is zero based on the fifth rule in (33). Using this result, we
have:

$\displaystyle\Pi^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\Sigma\left({\bm{\delta}}^{(t)}({\bm{x}}),{\bm{\delta}}^{(t^{\prime})}({\bm{x}})\right)$
$\displaystyle={\mathbb{E}}\left[[{\bm{\delta}}^{(t)}({\bm{x}})]_{i}\bm{\cdot}[{\bm{\delta}}^{(t^{\prime})}({\bm{x}}^{\prime})]_{i}\right]$
$\displaystyle=\sigma_{w}^{2}{\mathbb{E}}\left[[\phi^{\prime}({\bm{g}}^{(t)}({\bm{x}}))]_{i}\bm{\cdot}[{\bm{\delta}}^{(t+1)}({\bm{x}})]_{i}\bm{\cdot}[\phi^{\prime}({\bm{g}}^{(t^{\prime})}({\bm{x}}^{\prime}))]_{i}\bm{\cdot}[{\bm{\delta}}^{(t^{\prime}+1)}({\bm{x}}^{\prime})]_{i}\right]$
$\displaystyle=\sigma_{w}^{2}\mathop{\mathbb{E}}_{{\mathbf{z}}\sim\mathcal{N}(0,{\bm{K}}^{(t+1,t+1^{\prime})}({\bm{x}},{\bm{x}}^{\prime}))}\left[\phi^{\prime}({\textnormal{z}}_{1})\bm{\cdot}\phi^{\prime}({\textnormal{z}}_{2})\right]\bm{\cdot}{\mathbb{E}}\left[[{\bm{\delta}}^{(t+1)}({\bm{x}})]_{i}\bm{\cdot}[{\bm{\delta}}^{(t^{\prime}+1)}({\bm{x}}^{\prime})]_{i}\right]$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(t+1,t^{\prime}+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}\Pi^{(t+1,t^{\prime}+1)}({\bm{x}},{\bm{x}}^{\prime}).$

If $T^{\prime}-t^{\prime}=T-t$, then the the formula will lead to

$\displaystyle\Pi^{(T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle={\mathbb{E}}\left[[{\bm{\delta}}^{(T)}({\bm{x}})]_{i},[{\bm{\delta}}^{(T^{\prime})}({\bm{x}}^{\prime})]_{i}\right]$
$\displaystyle=\sigma_{v}^{2}{\mathbb{E}}\left[[{\bm{v}}]_{i}\bm{\cdot}[\phi^{\prime}({\bm{g}}^{(T)}({\bm{x}}))]_{i}\bm{\cdot}[{\bm{v}}]_{i}\bm{\cdot}[\phi^{\prime}({\bm{g}}^{(T^{\prime})}({\bm{x}}^{\prime}))]_{i}\right]$
$\displaystyle={\mathbb{E}}\left[[\phi^{\prime}({\bm{g}}^{(T)}({\bm{x}}))]_{i}\bm{\cdot}[\phi^{\prime}({\bm{g}}^{(T^{\prime})}({\bm{x}}^{\prime}))]_{i}\right]\bm{\cdot}{\mathbb{E}}\left[\left[{\bm{v}}\right]_{i}\left[{\bm{v}}\right]_{i}\right]$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(T+1,T+\tau+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}.$

Otherwise it will end to either of two cases for some $t^{\prime\prime}<T$ or
$T^{\prime}$ and by the fifth rule in (33) we have:

$\displaystyle\Sigma\left({\bm{\delta}}^{(t^{\prime\prime})}({\bm{x}}),{\bm{\delta}}^{(T^{\prime})}({\bm{x}})\right)=\Sigma\left(\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}^{\top}\Big{(}\phi^{\prime}({\bm{g}}^{(t^{\prime\prime})}({\bm{x}}))\odot{\bm{\delta}}^{(t^{\prime\prime}+1)}({\bm{x}}^{\prime})\Big{)},{\bm{v}}\odot\phi^{\prime}({\bm{g}}^{(T^{\prime})}({\bm{x}}))\right)=0$
$\displaystyle\Sigma\left({\bm{\delta}}^{(T)}({\bm{x}}),{\bm{\delta}}^{(t^{\prime\prime})}({\bm{x}})\right)=\Sigma\left({\bm{v}}\odot\phi^{\prime}({\bm{g}}^{(T)}({\bm{x}})),\frac{\sigma_{w}}{\sqrt{n}}{\mathbf{W}}^{\top}\Big{(}\phi^{\prime}({\bm{g}}^{(t^{\prime\prime})}({\bm{x}}^{\prime}))\odot{\bm{\delta}}^{(t^{\prime\prime}+1)}({\bm{x}}^{\prime})\Big{)}\right)=0.$

Without loss of generality, from now on assume $T^{\prime}<T$ and
$T^{\prime}-T=\tau$, the final formula for computing the backward gradients
becomes:

$\displaystyle\Pi^{(T,T+\tau)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(T+1,T+\tau+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}$
$\displaystyle\Pi^{(t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(t+1,t+\tau+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}\Pi^{(t+1,t+1+\tau)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle t\in[T-1]$
$\displaystyle\Pi^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=0$ $\displaystyle t^{\prime}-t\neq\tau$

Now we have derived the single layer RNTK. Recall that
$\theta=\mathrm{Vect}\big{[}\\{{\mathbf{W}},{\mathbf{U}},{\mathbf{b}},{\bm{v}}\\}\big{]}$
contains all of the network’s learnable parameters. As a result, we have:

$\displaystyle\nabla_{\theta}f_{\theta}({\bm{x}})=\mathrm{Vect}\big{[}\\{\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{W}}},\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}},\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{b}}},\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\bm{v}}}\\}\big{]}.$

As a result

$\displaystyle\left\langle\nabla_{\theta}f_{\theta}({\bm{x}}),\nabla_{\theta}f_{\theta}({\bm{x}}^{\prime})\right\rangle$
$\displaystyle=\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{W}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{W}}}\right\rangle+\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{U}}}\right\rangle+\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{b}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{b}}}\right\rangle$
$\displaystyle+\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\bm{v}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\bm{v}}}\right\rangle$

Where the gradients of output with respect to weights can be formulated as the
following compact form:

$\displaystyle\frac{\partial f_{\theta}({\bm{x}})}{\partial{\mathbf{W}}}$
$\displaystyle=\sum_{t=1}^{T}\left(\frac{1}{\sqrt{n}}{\bm{\delta}}^{(t)}({\bm{x}})\right)\bm{\cdot}\left(\frac{\sigma_{w}}{\sqrt{n}}{\bm{h}}^{(t-1)}({\bm{x}})\right)^{\top}$
$\displaystyle\frac{\partial f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}}$
$\displaystyle=\sum_{t=1}^{T}\left(\frac{1}{\sqrt{n}}{\bm{\delta}}^{(t)}({\bm{x}})\right)\bm{\cdot}\left(\frac{\sigma_{u}}{\sqrt{m}}{\bm{x}}_{t}\right)^{\top}$
$\displaystyle\frac{\partial f_{\theta}({\bm{x}})}{\partial{\mathbf{b}}}$
$\displaystyle=\sum_{t=1}^{T}\left(\frac{\sigma_{b}}{\sqrt{n}}{\bm{\delta}}^{(t)}({\bm{x}})\right)$
$\displaystyle\frac{\partial f_{\theta}({\bm{x}})}{\partial{\bm{v}}}$
$\displaystyle=\frac{\sigma_{v}}{\sqrt{n}}{\bm{h}}^{(T)}({\bm{x}}).$

As a result we have:

$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{W}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{W}}}\right\rangle$
$\displaystyle=\sum_{t^{\prime}=1}^{T^{\prime}}\sum_{t=1}^{T}\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(t)}({\bm{x}}),{\bm{\delta}}^{(t^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}\left(\frac{\sigma_{w}^{2}}{n}\left\langle{\bm{h}}^{(t-1)}({\bm{x}}),{\bm{h}}^{(t^{\prime}-1)}({\bm{x}}^{\prime})\right\rangle\right)$
$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{U}}}\right\rangle$
$\displaystyle=\sum_{t^{\prime}=1}^{T^{\prime}}\sum_{t=1}^{T}\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(t)}({\bm{x}}),{\bm{\delta}}^{(t^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}\left(\frac{\sigma_{u}^{2}}{m}\left\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{t^{\prime}}\right\rangle\right)$
$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{b}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{b}}}\right\rangle$
$\displaystyle=\sum_{t^{\prime}=1}^{T^{\prime}}\sum_{t=1}^{T}\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(t)}({\bm{x}}),{\bm{\delta}}^{(t^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}\sigma_{b}^{2}$
$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\bm{v}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\bm{v}}}\right\rangle$
$\displaystyle=\left(\frac{\sigma_{v}^{2}}{n}\left\langle{\bm{h}}^{(T)}({\bm{x}}),{\bm{h}}^{(T^{\prime})}({\bm{x}}^{\prime})\right\rangle\right).$

Remember that for any two $G$-var
${\mathbb{E}}\left[[{\bm{g}}]_{i}[{\bm{g}}^{\prime}]_{i}\right]$ is
independent of index $i$. Therefore,

$\displaystyle\frac{1}{n}\left\langle{\bm{h}}^{(t-1)}({\bm{x}}),{\bm{h}}^{(t^{\prime}-1)}({\bm{x}}^{\prime})\right\rangle$
$\displaystyle\rightarrow\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\big{]}\hskip
20.0ptt>1$
$\displaystyle\frac{1}{n}\left\langle{\bm{h}}^{(0)}({\bm{x}}),{\bm{h}}^{(0)}({\bm{x}}^{\prime})\right\rangle$
$\displaystyle\rightarrow\sigma_{h}^{2}.$

Hence, by summing the above terms in the infinite-width limit we get

$\displaystyle\left\langle\nabla_{\theta}f_{\theta}({\bm{x}}),\nabla_{\theta}f_{\theta}({\bm{x}}^{\prime})\right\rangle\rightarrow\left(\sum_{t^{\prime}=1}^{T^{\prime}}\sum_{t=1}^{T}\Pi^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Sigma^{(t,t^{\prime})}({\bm{x}}^{\prime},{\bm{x}}^{\prime})\right)+\mathcal{K}({\bm{x}},{\bm{x}}^{\prime}).$
(37)

Since $\Pi^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})=0$ for
$t^{\prime}-t\neq\tau$ it is simplified to

$\displaystyle\left\langle\nabla_{\theta}f_{\theta}({\bm{x}}),\nabla_{\theta}f_{\theta}({\bm{x}}^{\prime})\right\rangle=\left(\sum_{t=1}^{T}\Pi^{(t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Sigma^{(t,t+\tau)}({\bm{x}}^{\prime},{\bm{x}}^{\prime})\right)+\mathcal{K}({\bm{x}},{\bm{x}}^{\prime}).$

Multi-dimensional output. For $f_{\theta}({\bm{x}})\in\mathbb{R}^{d}$, the
$i$-th output for $i\in[d]$ is obtained via

$\displaystyle\left[f_{\theta}({\bm{x}})\right]_{i}=\frac{\sigma_{v}}{\sqrt{n}}{\bm{v}}^{\top}_{i}{\bm{h}}^{(T)}({\bm{x}}),$

where ${\bm{v}}_{i}$ is independent of ${\bm{v}}_{j}$ for $i\neq j$. As a
result, for The RNTK $\Theta({\bm{x}},{\bm{x}}^{\prime})\in\mathbb{R}^{d\times
d}$ for multi-dimensional output we have

$\displaystyle[\Theta({\bm{x}},{\bm{x}}^{\prime})]_{i,j}$
$\displaystyle=\left\langle\nabla_{\theta}\left[f_{\theta}({\bm{x}})\right]_{i},\nabla_{\theta}\left[f_{\theta}({\bm{x}}^{\prime})\right]_{j}\right\rangle$

For $i=j$, the kernel is the same as computed in (37) and we denote it as

$\displaystyle\left\langle\nabla_{\theta}\left[f_{\theta}({\bm{x}})\right]_{i},\nabla_{\theta}\left[f_{\theta}({\bm{x}}^{\prime})\right]_{i}\right\rangle=\Theta^{(T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime}).$

For $i\neq j$, since ${\bm{v}}_{i}$ is independent of ${\bm{v}}_{j}$,
$\Pi^{(T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$ and all the backward pass
gradients become zero, so

$\displaystyle\left\langle\nabla_{\theta}\left[f_{\theta}({\bm{x}})\right]_{i},\nabla_{\theta}\left[f_{\theta}({\bm{x}}^{\prime})\right]_{j}\right\rangle=0\hskip
20.0pti\neq j$

which gives us the following formula

$\displaystyle\Theta({\bm{x}},{\bm{x}}^{\prime})=\Theta^{(T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\otimes{\bm{I}}_{d}.$

This concludes the proof for Theorem 1 for single-layer case.

###  B.3 Proof for Theorem 1: Multi-Layer Case

Now we drive the RNTK for multi-layer RNTK. We will only study single output
case and the generalization to multi-dimensional case is identical as the
single layer case. The set of equations for calculation of the output of a
$L$-layer RNN for ${\bm{x}}=\\{{\bm{x}}_{t}\\}_{t=1}^{T}$ are

$\displaystyle{\bm{g}}^{\left(\ell,t\right)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{w}^{\ell}}{\sqrt{n}}{\mathbf{W}}^{\left(\ell\right)}{\bm{h}}^{\left(\ell,t-1\right)}({\bm{x}})+\frac{\sigma_{u}^{\ell}}{\sqrt{m}}{\mathbf{U}}^{\left(\ell\right)}{\bm{x}}_{t}+\sigma_{b}^{\ell}{\mathbf{b}}^{\left(\ell\right)}\hskip
10.0pt$ $\displaystyle\ell=1$
$\displaystyle{\bm{g}}^{\left(\ell,t\right)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{w}^{\ell}}{\sqrt{n}}{\mathbf{W}}^{\left(\ell\right)}{\bm{h}}^{\left(\ell,t-1\right)}({\bm{x}})+\frac{\sigma_{u}^{\ell}}{\sqrt{n}}{\mathbf{U}}^{\left(\ell\right)}{\bm{h}}^{\left(\ell-1,t\right)}({\bm{x}})+\sigma_{b}^{\ell}{\mathbf{b}}^{\left(\ell\right)}\hskip
10.0pt$ $\displaystyle\ell>1$
$\displaystyle{\bm{h}}^{\left(\ell,t\right)}({\bm{x}})$
$\displaystyle=\phi\left({\bm{g}}^{\left(\ell,t\right)}({\bm{x}})\right)$
$\displaystyle f_{\theta}({\bm{x}})$
$\displaystyle=\frac{\sigma_{v}}{\sqrt{n}}{\bm{v}}^{\top}{\bm{h}}^{(L,T)}({\bm{x}})$

The forward pass kernels for the first layer is the same as calculated in B.2.
For $\ell\geq 2$ we have:

$\displaystyle\Sigma^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\Sigma({\bm{g}}^{\left(\ell,t\right)}({\bm{x}}),{\bm{g}}^{(\ell,t^{\prime})}({\bm{x}}^{\prime}))$
$\displaystyle=\Sigma\left(\frac{\sigma_{w}^{\ell}}{\sqrt{n}}{\mathbf{W}}^{(\ell)}{\bm{h}}^{(\ell,t-1)}({\bm{x}}),\frac{\sigma_{w}^{\ell}}{\sqrt{n}}{\mathbf{W}}^{(\ell)}{\bm{h}}^{(\ell,t^{\prime}-1)}({\bm{x}}^{\prime})\right)$
$\displaystyle+\Sigma\left(\frac{\sigma_{u}^{\ell}}{\sqrt{n}}{\mathbf{U}}^{(\ell)}{\bm{h}}^{\left(\ell-1,t\right)}({\bm{x}}),\frac{\sigma_{u}^{\ell}}{\sqrt{n}}{\mathbf{U}}^{(\ell)}{\bm{h}}^{\left(\ell-1,t^{\prime}\right)}({\bm{x}}^{\prime})\right)+\Sigma^{\rm
in}\left(\sigma_{b}^{\ell}{\bm{b}}^{(\ell)},\sigma_{b}^{\ell}{\bm{b}}^{(\ell)}\right)$
$\displaystyle=(\sigma_{w}^{\ell})^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\big{]}+(\sigma_{u}^{\ell})^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(\ell-1,t+1,t^{\prime}+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}+(\sigma_{b}^{\ell})^{2},$

where

$\displaystyle{\bm{K}}^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\left[{\begin{array}[]{cc}{\bm{\Sigma}}^{(\ell,t-1,t-1)}({\bm{x}},{\bm{x}})&{\bm{\Sigma}}^{(\ell,t-1,t^{\prime}-1)}({\bm{x}},{\bm{x}}^{\prime})\\\
{\bm{\Sigma}}^{(\ell,t-1,t^{\prime}-1)}({\bm{x}},{\bm{x}}^{\prime})&{\bm{\Sigma}}^{(\ell,t^{\prime}-1,t^{\prime}-1)}({\bm{x}}^{\prime},{\bm{x}}^{\prime})\\\
\end{array}}\right]\,,$

and $\Sigma^{\rm in}$ is defined in (B.2). For the first first time step we
have:

$\displaystyle\Sigma^{(\ell,1,1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=(\sigma_{w}^{\ell})^{2}\sigma_{h}^{2}1_{({\bm{x}}={\bm{x}}^{\prime})}+(\sigma_{u}^{\ell})^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(\ell,2,2)}({\bm{x}},{\bm{x}}^{\prime})\big{]}+(\sigma_{b}^{\ell})^{2}\,,$
$\displaystyle\Sigma^{(\ell,t,1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=(\sigma_{u}^{\ell})^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(\ell,t+1,2)}({\bm{x}},{\bm{x}}^{\prime})\big{]}+(\sigma_{b}^{\ell})^{2}\,,$
$\displaystyle\Sigma^{(\ell,1,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=(\sigma_{u}^{\ell})^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(\ell,2,t^{\prime}+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}+(\sigma_{b}^{\ell})^{2}\,.$

And the output layer

$\displaystyle\mathcal{K}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(L,T+1,T^{\prime}+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}.$

Note that because of using new weights at each layer we get

$\displaystyle\Sigma({\bm{g}}^{(\ell,t)}({\bm{x}}),{\bm{g}}^{(\ell^{\prime},t^{\prime})})({\bm{x}}))=0\hskip
20.0pt\ell\neq\ell^{\prime}$

Now we calculate the backward pass kernels in multi-layer RNTK. The gradients
at the last layer is calculated via

$\displaystyle{\bm{\delta}}^{(L,T)}({\bm{x}})$
$\displaystyle=\sigma_{v}{\bm{v}}\odot\phi^{\prime}({\bm{g}}^{(L,T)}({\bm{x}})).$

In the last hidden layer for different time steps we have

$\displaystyle{\bm{\delta}}^{(L,t)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{w}^{L}}{\sqrt{n}}\left({\mathbf{W}}^{(L)}\right)^{\top}\Big{(}\phi^{\prime}({\bm{g}}^{(L,t)}({\bm{x}}))\odot{\bm{\delta}}^{(L,t+1)}({\bm{x}})\Big{)}\hskip
20.0ptt\in[T-1]$

In the last time step for different hidden layers we have

$\displaystyle{\bm{\delta}}^{(\ell,T)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{u}^{\ell+1}}{\sqrt{n}}\left({\mathbf{U}}^{(\ell+1)}\right)^{\top}\Big{(}\phi^{\prime}({\bm{g}}^{(\ell,T)}({\bm{x}}))\odot{\bm{\delta}}^{(\ell+1,T)}({\bm{x}})\Big{)}\hskip
20.0pt\ell\in[L-1]$

At the end for the other layers we have

$\displaystyle{\bm{\delta}}^{(\ell,t)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{w}^{\ell}}{\sqrt{n}}\left({\mathbf{W}}^{(\ell)}\right)^{\top}\Big{(}\phi^{\prime}({\bm{g}}^{(\ell,t)}({\bm{x}}))\odot{\bm{\delta}}^{(\ell,t+1)}({\bm{x}})\Big{)}$
$\displaystyle+\frac{\sigma_{u}^{\ell+1}}{\sqrt{n}}\left({\mathbf{U}}^{(\ell+1)}\right)^{\top}\Big{(}\phi^{\prime}({\bm{g}}^{(\ell,t)}({\bm{x}}))\odot{\bm{\delta}}^{(\ell+1,t)}({\bm{x}})\Big{)}\hskip
20.0pt\ell\in[L-1],t\in[T-1]$

The recursive formula for the
$\Pi^{(L,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$ is the same as the single
layer, and it is non-zero for $t^{\prime}-t=T^{\prime}-T=\tau$. As a result we
have

$\displaystyle\Pi^{(L,T,T+\tau)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(L,T+1,T+\tau+1)}\big{]}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle\Pi^{(L,t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=(\sigma_{w}^{L})^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(L,t+1,t+\tau+1)}\big{]}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Pi^{(L,t+1,t+\tau+1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle t\in[T-1]$
$\displaystyle\Pi^{(L,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=0$ $\displaystyle t^{\prime}-t\neq\tau$ (38)

Similarly by using the same course of arguments used in the single layer
setting, for the last time step we have

$\displaystyle\Pi^{(\ell,T,T+\tau)}({\bm{x}},{\bm{x}}^{\prime})=(\sigma^{\ell+1}_{u})^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(\ell,T+1,T+\tau+1)}\big{]}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Pi^{(\ell+1,T,T+\tau)}({\bm{x}},{\bm{x}}^{\prime})\hskip
15.0pt\ell\in[L-1]$

For the other layers we have

$\displaystyle\Pi^{(\ell,t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=(\sigma_{w}^{\ell})^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(\ell,t+1,t+\tau+1)}\big{]}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Pi^{(\ell,t+1,t+\tau+1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle+(\sigma^{\ell+1}_{u})^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(\ell,t+1,t+\tau+1)}\big{]}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Pi^{(\ell+1,t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime}).$

For $t^{\prime}-t\neq\tau$ the recursion continues until it reaches
$\Pi^{(L,T,t^{\prime\prime})}({\bm{x}},{\bm{x}}^{\prime}),t^{\prime\prime}<T^{\prime}$
or
$\Pi^{(L,t^{\prime\prime},T^{\prime})}({\bm{x}},{\bm{x}}^{\prime}),t^{\prime\prime}<T$
and as a result based on (38) we get

$\displaystyle\Pi^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})=0\hskip
20.0ptt^{\prime}-t\neq\tau$ (39)

For $t^{\prime}-t=\tau$ it leads to
$\Pi^{(L,T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$ and has a non-zero value.
Now we derive RNTK for multi-layer:

$\displaystyle\left\langle\nabla_{\theta}f_{\theta}({\bm{x}}),\nabla_{\theta}f_{\theta}({\bm{x}}^{\prime})\right\rangle$
$\displaystyle=\sum_{\ell=1}^{L}\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{W}}^{(\ell)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{W}}^{(\ell)}}\right\rangle+\sum_{\ell=1}^{L}\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}^{(\ell)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{U}}^{(\ell)}}\right\rangle$
$\displaystyle+\sum_{\ell=1}^{L}\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{b}}^{(\ell)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{b}}^{(\ell)}}\right\rangle+\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\bm{v}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\bm{v}}}\right\rangle,$

where

$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{W}}^{(\ell)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{W}}^{(\ell)}}\right\rangle$
$\displaystyle=\sum_{t^{\prime}=1}^{T^{\prime}}\sum_{t=1}^{T}\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(\ell,t)}({\bm{x}}),{\bm{\delta}}^{(\ell,t^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}\left(\frac{(\sigma_{w}^{\ell})^{2}}{n}\left\langle{\bm{h}}^{(\ell,t-1)}({\bm{x}}),{\bm{h}}^{(\ell,t^{\prime}-1)}({\bm{x}}^{\prime})\right\rangle\right)$
$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}^{(\ell)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{U}}^{(\ell)}}\right\rangle$
$\displaystyle=\sum_{t^{\prime}=1}^{T^{\prime}}\sum_{t=1}^{T}\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(\ell,t)}({\bm{x}}),{\bm{\delta}}^{(\ell,t^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}\left(\frac{(\sigma_{u}^{\ell})^{2}}{m}\left\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{t^{\prime}}\right\rangle\right)\hskip
20.0pt\ell=1$ $\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}^{(\ell)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{U}}^{(\ell)}}\right\rangle$
$\displaystyle=\sum_{t^{\prime}=1}^{T^{\prime}}\sum_{t=1}^{T}\bigg{[}\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(\ell,t)}({\bm{x}}),{\bm{\delta}}^{(\ell,t^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)$
$\displaystyle\hskip
50.0pt\bm{\cdot}\left(\frac{(\sigma_{u}^{\ell})^{2}}{n}\left\langle{\bm{h}}^{(\ell-1,t)}({\bm{x}}),{\bm{h}}^{(\ell-1,t^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)\bigg{]}\hskip
47.0pt\ell>1$ $\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{b}}^{(\ell)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{b}}^{(\ell)}}\right\rangle$
$\displaystyle=\sum_{t^{\prime}=1}^{T^{\prime}}\sum_{t=1}^{T}\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(\ell,t)}({\bm{x}}),{\bm{\delta}}^{(\ell,t^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}(\sigma_{b}^{\ell})^{2}$
$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\bm{v}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\bm{v}}}\right\rangle$
$\displaystyle=\left(\frac{\sigma_{v}^{2}}{n}\left\langle{\bm{h}}^{(T)}({\bm{x}}),{\bm{h}}^{(T^{\prime})}({\bm{x}}^{\prime})\right\rangle\right)$

Summing up all the terms and replacing the inner product of vectors with their
expectations we get

$\displaystyle\left\langle\nabla_{\theta}f_{\theta}({\bm{x}}),\nabla_{\theta}f_{\theta}({\bm{x}}^{\prime})\right\rangle=\Theta^{(L,T,T^{\prime})}=\left(\sum_{\ell=1}^{L}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T^{\prime}}\Pi^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Sigma^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\right)+\mathcal{K}({\bm{x}},{\bm{x}}^{\prime}).$

By (39), we can simplify to

$\displaystyle\Theta^{(L,T,T^{\prime})}=\left(\sum_{\ell=1}^{L}\sum_{t=1}^{T}\Pi^{(\ell,t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Sigma^{(\ell,t,t+\tau)}({\bm{x}},{\bm{x}}^{\prime})\right)+\mathcal{K}({\bm{x}},{\bm{x}}^{\prime}).$

For multi-dimensional output it becomes

$\displaystyle\Theta({\bm{x}},{\bm{x}}^{\prime})=\Theta^{(L,T,T^{\prime})}({\bm{x}},{\bm{x}}^{\prime})\otimes{\bm{I}}_{d}.$

This concludes the proof for Theorem 1 for the multi-layer case.

###  B.4 Proof for Theorem 3: Weight-Untied RNTK

The architecture of a weight-untied single layer RNN is

$\displaystyle{\bm{g}}^{\left(t\right)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{w}}{\sqrt{m}}{\mathbf{W}}^{(t)}{\bm{h}}^{\left(t-1\right)}({\bm{x}})+\frac{\sigma_{u}}{\sqrt{n}}{\mathbf{U}}^{(t)}{\bm{x}}_{t}+\sigma_{b}{\mathbf{b}}^{(t)}$
$\displaystyle{\bm{h}}^{\left(t\right)}({\bm{x}})$
$\displaystyle=\phi\left({\bm{g}}^{\left(t\right)}({\bm{x}})\right)$
$\displaystyle f_{\theta}({\bm{x}})$
$\displaystyle=\frac{\sigma_{v}}{\sqrt{n}}{\bm{v}}^{\top}{\bm{h}}^{(T)}({\bm{x}})$

Where we use new weights at each time step and we index it by time. Like
previous sections, we first derive the forward pass kernels for two same
length data
${\bm{x}}=\\{{\bm{x}}_{t}\\}_{t=1}^{T}$,${\bm{x}}=\\{{\bm{x}}^{\prime}_{t^{\prime}}\\}_{t^{\prime}=1}^{T}$

$\displaystyle\Sigma^{(t,t)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi}\big{[}{\bm{K}}^{(t,t)}({\bm{x}},{\bm{x}}^{\prime})\big{]}+\frac{\sigma_{u}^{2}}{m}\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{t}\rangle+\sigma_{b}^{2}.$
$\displaystyle\Sigma^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=0\hskip 20.0ptt\neq t^{\prime}$

Since we are using same weight at the same time step,
$\Sigma^{(t,t)}({\bm{x}},{\bm{x}}^{\prime})$ can be written as a function of
the previous kernel, which is exactly as the weight-tied RNN. However for
different length, it becomes zero as a consequence of using different weights,
unlike weight-tied which has non-zero value. The kernel of the first time step
and output is also the same as weight-tied RNN. For the gradients we have:

$\displaystyle{\bm{\delta}}^{(T)}({\bm{x}})$
$\displaystyle=\sigma_{v}{\bm{v}}\odot\phi^{\prime}({\bm{g}}^{(T)}({\bm{x}}))$
$\displaystyle{\bm{\delta}}^{(t)}({\bm{x}})$
$\displaystyle=\frac{\sigma_{w}}{\sqrt{n}}({\mathbf{W}}^{(t+1)})^{\top}\Big{(}\phi^{\prime}({\bm{g}}^{(t)}({\bm{x}}))\odot{\bm{\delta}}^{(t+1)}({\bm{x}})\Big{)}\hskip
20.0ptt\in[T-1]$

For $t^{\prime}=t$ we have:

$\displaystyle\Pi^{(t,t)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{w}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(t+1,t+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}\Pi^{(t+1,t+1)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle\Pi^{(t,t)}({\bm{x}},{\bm{x}}^{\prime})$
$\displaystyle=\sigma_{v}^{2}\mathrm{V}_{\phi^{\prime}}\big{[}{\bm{K}}^{(T+1,T+\tau+1)}({\bm{x}},{\bm{x}}^{\prime})\big{]}.$

Due to using different weights for $t\neq t^{\prime}$, we can immediately
conclude that $\Pi^{(t,t^{\prime})}({\bm{x}},{\bm{x}}^{\prime})=0$. This set
of calculation is exactly the same as the weight-tied case when $\tau=T-T=0$.
Finally, with
$\theta=\mathrm{Vect}\big{[}\\{\\{{\mathbf{W}}^{(t)},{\mathbf{U}}^{(t)},{\mathbf{b}}^{(t)}\\}_{t=1}^{T},{\bm{v}}\\}\big{]}$
we have

$\displaystyle\left\langle\nabla_{\theta}f_{\theta}({\bm{x}}),\nabla_{\theta}f_{\theta}({\bm{x}}^{\prime})\right\rangle$
$\displaystyle=\sum_{t=1}^{T}\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{W}}^{(t)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{W}}^{(t)}}\right\rangle+\sum_{t=1}^{T}\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}^{(t)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{U}}^{(t)}}\right\rangle$
$\displaystyle+\sum_{t=1}^{T}\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{b}}^{(t)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{b}}^{(t)}}\right\rangle+\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\bm{v}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\bm{v}}},\right\rangle$

with

$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{W}}^{(t)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{W}}^{(t)}}\right\rangle$
$\displaystyle=\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(t)}({\bm{x}}),{\bm{\delta}}^{(t)}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}\left(\frac{\sigma_{w}^{2}}{n}\left\langle{\bm{h}}^{(t-1)}({\bm{x}}),{\bm{h}}^{(t-1)}({\bm{x}}^{\prime})\right\rangle\right)$
$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{U}}^{(t)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{U}}^{(t)}}\right\rangle$
$\displaystyle=\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(t)}({\bm{x}}),{\bm{\delta}}^{(t)}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}\left(\frac{\sigma_{u}^{2}}{m}\left\langle{\bm{x}}_{t},{\bm{x}}^{\prime}_{t}\right\rangle\right)$
$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\mathbf{b}}^{(t)}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\mathbf{b}}^{(t)}}\right\rangle$
$\displaystyle=\left(\frac{1}{n}\left\langle{\bm{\delta}}^{(t)}({\bm{x}}),{\bm{\delta}}^{(t)}({\bm{x}}^{\prime})\right\rangle\right)\bm{\cdot}\sigma_{b}^{2}$
$\displaystyle\left\langle\frac{\partial
f_{\theta}({\bm{x}})}{\partial{\bm{v}}},\frac{\partial
f_{\theta}({\bm{x}}^{\prime})}{\partial{\bm{v}}}\right\rangle$
$\displaystyle=\left(\frac{\sigma_{v}^{2}}{n}\left\langle{\bm{h}}^{(T)}({\bm{x}}),{\bm{h}}^{(T^{\prime})}({\bm{x}}^{\prime})\right\rangle\right).$

As a result we obtain

$\displaystyle\left\langle\nabla_{\theta}f_{\theta}({\bm{x}}),\nabla_{\theta}f_{\theta}({\bm{x}}^{\prime})\right\rangle=\left(\sum_{t=1}^{T}\Pi^{(t,t)}({\bm{x}},{\bm{x}}^{\prime})\bm{\cdot}\Sigma^{(t,t)}({\bm{x}}^{\prime},{\bm{x}}^{\prime})\right)+\mathcal{K}({\bm{x}},{\bm{x}}^{\prime}),$

same as the weight-tied RNN when $\tau=0$. This concludes the proof for
Theorem 3.

###  B.5 Analytical Formula for $\mathrm{V}_{\phi}[K]$

For any positive definite matrix
${\bm{K}}=\left[{\begin{array}[]{cc}K_{1}&K_{3}\\\ K_{3}&K_{2}\\\
\end{array}}\right]$ we have:

  * •

$\phi=\rm ReLU$ [12]

$\displaystyle\mathrm{V}_{\phi}[{\bm{K}}]$
$\displaystyle=\frac{1}{2\pi}\left(c(\pi-\rm
arccos(c))+\sqrt{1-c^{2}})\right)\sqrt{K_{1}K_{2}},$
$\displaystyle\mathrm{V}_{\phi^{\prime}}[{\bm{K}}]$
$\displaystyle=\frac{1}{2\pi}(\pi-\rm arccos(c)).$

where $c=K_{3}/\sqrt{K_{1}K_{2}}$

  * •

$\phi=\rm erf$ [30]

$\displaystyle\mathrm{V}_{\phi}[{\bm{K}}]$ $\displaystyle=\frac{2}{\pi}\rm
arcsin\left(\frac{2K_{3}}{\sqrt{(1+2K_{1})(1+2K_{3})}}\right),$
$\displaystyle\mathrm{V}_{\phi^{\prime}}[{\bm{K}}]$
$\displaystyle=\frac{4}{\pi\sqrt{(1+2K_{1})(1+2K_{2})-4K_{3}^{2}}}.$

##  Appendix C Proof for Theorem 2: RNTK Convergence after Training

To prove theorem 2, we use the strategy used in [29] which relies on the the
local lipschitzness of the network Jacobian
${\bm{J}}(\theta,\mathcal{X})=\nabla_{\theta}f_{\theta}({\bm{x}})\in\mathbb{R}^{|\mathcal{X}|d\times|\theta|}$
at initialization.

###### Definition 1

The Jacobian of a neural network is local lipschitz at NTK initialization
$\left(\theta_{0}\sim\mathcal{N}(0,1)\right)$ if there is constant $K>0$ for
every $C$ such that

$\displaystyle\begin{dcases}\|{\bm{J}}(\theta,\mathcal{X})\|_{F}<K\\\
\|{\bm{J}}(\theta,\mathcal{X})-{\bm{J}}(\tilde{\theta},\mathcal{X})\|_{F}<K\|\theta-\tilde{\theta}\|\end{dcases},\hskip
30.0pt\forall\,\,\theta,\tilde{\theta}\in B(\theta_{0},R)$

where

$\displaystyle B(\theta,R):=\\{\theta:\|\theta_{0}-\theta\|<R\\}.$

###### Theorem 4

Assume that the network Jacobian is local lipschitz with high probability and
the empirical NTK of the network converges in probability at initialization
and it is positive definite over the input set. For $\epsilon>0$, there exists
$N$ such that for $n>N$ when applying gradient flow with
$\eta<2\left(\lambda_{\mathrm{min}}(\Theta(\mathcal{X},\mathcal{X}))+\lambda_{\mathrm{max}}(\Theta(\mathcal{X},\mathcal{X})\right)^{-1}$
with probability at least $(1-\epsilon)$ we have:

$\displaystyle\underset{s}{\mathrm{sup}}\frac{\|\theta_{s}-\theta_{0}\|_{2}}{\sqrt{n}},\underset{s}{\mathrm{sup}}\|\widehat{\Theta}_{s}(\mathcal{X},\mathcal{X})-\widehat{\Theta}_{0}(\mathcal{X},\mathcal{X})\|=\mathcal{O}\left(\frac{1}{\sqrt{n}}\right).$

Proof: See [29]

Theorem 4 holds for any network architecture and any cost function and it was
used in [29] to show the stability of NTK for MLP during training.

Here we extend the results for RNTK by proving that the Jacobian of a multi-
layer RNN under NTK initialization is local lipschitz with high probability.

To prove it, first, we prove that for any two points $\theta,\tilde{\theta}\in
B(\theta_{0},R)$ there exists constant $K_{1}$ such that

$\displaystyle\|{\bm{g}}^{(\ell,t)}({\bm{x}})\|_{2},\|\delta^{(\ell,t)}({\bm{x}})\|_{2}\leq
K_{1}\sqrt{n}$ (40)
$\displaystyle\|{\bm{g}}^{(\ell,t)}({\bm{x}})-\tilde{{\bm{g}}}^{(\ell,t)}({\bm{x}})\|_{2},\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})-\tilde{{\bm{\delta}}}^{(\ell,t)}({\bm{x}})\|_{2}\leq\|\bar{\theta}-\tilde{\theta}\|\leq
K_{1}\sqrt{n}\|\theta-\tilde{\theta}\|.$ (41)

To prove (40) and (41) we use the following lemmas.333See
math.uci.edu/~rvershyn/papers/HDP-book/HDP-book.pdf for proofs

###### Lemma 1

Let $A\in\mathbb{R}^{n\times m}$ be a random matrix whose entries are
independent standard normal random variables. Then for every $t\geq 0$, with
probability at least $1-e^{(-ct^{2})}$ for some constant $c$ we have:

$\displaystyle\|A\|_{2}\leq\sqrt{m}+\sqrt{n}+t.$

###### Lemma 2

Let $a\in\mathbb{R}^{n}$ be a random vector whose entries are independent
standard normal random variables. Then for every $t\geq 0$, with probability
at least $1-e^{(-ct^{2})}$ for some constant $c$ we have:

$\displaystyle\|a\|_{2}\leq\sqrt{n}+\sqrt{t}.$

Setting $t=\sqrt{n}$ for any $\theta\in R(\theta_{0},R)$. With high
probability, we get:

$\displaystyle\|{\mathbf{W}}^{(\ell)}\|_{2},\|{\mathbf{U}}^{(\ell)}\|_{2}\leq
3\sqrt{n},\hskip 10.0pt\|{\mathbf{b}}^{\ell}\|_{2}\leq 2\sqrt{n},\hskip
10.0pt\|{\mathbf{h}}^{(\ell,0)}({\bm{x}})\|_{2}\leq 2\sigma_{h}\sqrt{n}.$

We also assume that there exists some finite constant $C$ such that

$\displaystyle|\phi(x)|<C|x|,\hskip
10.0pt|\phi(x)-\phi(x^{\prime})|<C|x-x^{\prime}|,\hskip
10.0pt|\phi^{\prime}(x)|<C,\hskip
10.0pt,|\phi^{\prime}(x)-\phi^{\prime}(x^{\prime})|<C|x-x^{\prime}|.$

The proof is obtained by induction. From now on assume that all inequalities
in (40) and (41) holds with some $k$ for the previous layers. We have

$\displaystyle\|{\bm{g}}^{(\ell,t)}({\bm{x}})\|_{2}$
$\displaystyle=\|\frac{\sigma_{w}^{\ell}}{\sqrt{n}}{\mathbf{W}}^{\left(\ell\right)}{\bm{h}}^{\left(\ell,t-1\right)}({\bm{x}})+\frac{\sigma_{u}^{\ell}}{\sqrt{n}}{\mathbf{U}}^{\left(\ell\right)}\
{\bm{h}}^{\left(\ell-1,t\right)}({\bm{x}})+\sigma_{b}^{\ell}{\mathbf{b}}^{\left(\ell\right)}\|_{2}$
$\displaystyle\leq\frac{\sigma_{w}^{\ell}}{\sqrt{n}}\|{\mathbf{W}}^{\left(\ell\right)}\|_{2}\|\phi\left({\bm{g}}^{\left(\ell,t-1\right)}({\bm{x}})\right)\|_{2}+\frac{\sigma_{u}^{\ell}}{\sqrt{n}}\|{\mathbf{U}}^{\left(\ell\right)}\|_{2}\|\phi\left({\bm{g}}^{\left(\ell-1,t\right)}({\bm{x}})\right)\|_{2}+\sigma_{b}^{\ell}\|{\mathbf{b}}^{\left(\ell\right)}\|_{2}$
$\displaystyle\leq\left(3\sigma_{w}^{\ell}Ck+3\sigma_{u}^{\ell}Ck+2\sigma_{b}\right)\sqrt{n}.$

And the proof for (40) and (41) is completed by showing that the first layer
is bounded

$\displaystyle\|{\bm{g}}^{(1,1)}({\bm{x}})\|_{2}$
$\displaystyle=\|\frac{\sigma_{w}^{1}}{\sqrt{n}}{\mathbf{W}}^{\left(\ell\right)}{\bm{h}}^{\left(1,0\right)}({\bm{x}})+\frac{\sigma_{u}^{1}}{\sqrt{m}}{\mathbf{U}}^{\left(\ell\right)}\
{\bm{x}}_{1}+\sigma_{b}^{1}{\mathbf{b}}^{\left(1\right)}\|_{2}$
$\displaystyle\leq(3\sigma_{w}^{1}\sigma_{h}+\frac{3\sigma_{u}}{\sqrt{m}}\|{\bm{x}}_{1}\|_{2}+2\sigma_{b})\sqrt{n}.$

For the gradient of first layer we have

$\displaystyle\|{\bm{\delta}}^{(L,T)}({\bm{x}})\|_{2}$
$\displaystyle=\|\sigma_{v}{\bm{v}}\odot\phi^{\prime}({\bm{g}}^{(L,T)}({\bm{x}}))\|_{2}$
$\displaystyle\leq\sigma_{v}\|{\bm{v}}\|_{2}\|\phi^{\prime}({\bm{g}}^{(L,T)}({\bm{x}}))\|_{\infty}$
$\displaystyle=2\sigma_{v}C\sqrt{n}.$

And similarly we have

$\displaystyle\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})\|$
$\displaystyle\leq\left(3\sigma_{w}Ck^{\prime}+3\sigma_{u}Ck^{\prime}\right)\sqrt{n}.$

For $\theta,\tilde{\theta}\in B(\theta_{0},R)$ we have

$\displaystyle\|{\bm{g}}^{(1,1)}({\bm{x}})-\tilde{{\bm{g}}}^{(1,1)}({\bm{x}})\|_{2}$
$\displaystyle=\|\frac{\sigma_{w}^{1}}{\sqrt{n}}({\mathbf{W}}^{(1)}-\tilde{{\mathbf{W}}}^{(1)}){\bm{h}}^{(1,0)}({\bm{x}})+\frac{\sigma_{u}^{1}}{\sqrt{m}}({\mathbf{U}}^{(1)}-\tilde{{\mathbf{U}}}^{(1)}){\bm{h}}^{(1,0)}({\bm{x}})\|_{2}$
$\displaystyle\leq\left(3\sigma_{w}^{1}\sigma_{h}+\frac{3\sigma_{u}^{1}}{m}\|{\bm{x}}_{1}\|_{2}\right)\|\theta-\tilde{\theta}\|_{2}\sqrt{n}.$
$\displaystyle\|{\bm{g}}^{(\ell,t)}({\bm{x}})-\tilde{{\bm{g}}}^{(\ell,t)}({\bm{x}})\|_{2}$
$\displaystyle\leq\|\phi({\bm{g}}^{(\ell,t-1)}({\bm{x}}))\|_{2}\|\frac{\sigma_{w}^{\ell}}{\sqrt{n}}({\mathbf{W}}^{(\ell)}-\tilde{{\mathbf{W}}}^{(\ell)})\|_{2}$
$\displaystyle+\|\frac{\sigma_{w}^{\ell}}{\sqrt{n}}\tilde{{\mathbf{W}}}^{(\ell)}\|_{2}\|\phi({\bm{g}}^{(\ell,t-1)}({\bm{x}}))-\phi(\tilde{{\bm{g}}}^{(\ell,t-1)}({\bm{x}}))\|_{2}$
$\displaystyle+\|\phi({\bm{g}}^{(\ell-1,t)}({\bm{x}}))\|_{2}\|\frac{\sigma_{u}^{\ell}}{\sqrt{n}}({\mathbf{U}}^{(\ell)}-\tilde{{\mathbf{U}}}^{(\ell)})\|_{2}$
$\displaystyle+\|\frac{\sigma_{u}^{\ell}}{\sqrt{n}}\tilde{{\mathbf{U}}}^{(\ell)}\|_{2}\|\phi({\bm{g}}^{(\ell-1,t)}({\bm{x}}))-\phi(\tilde{{\bm{g}}}^{(\ell-1,t)}({\bm{x}}))\|_{2}+\sigma_{b}\|{\mathbf{b}}^{(\ell)}-\tilde{{\mathbf{b}}}^{(\ell)}\|$
$\displaystyle\leq(k\sigma_{w}^{\ell}+3\sigma_{w}^{\ell}Ck+k\sigma_{u}^{\ell}+3\sigma_{u}^{\ell}Ck+\sigma_{b})\|\theta-\tilde{\theta}\|_{2}\sqrt{n}.$

For gradients we have

$\displaystyle\|{\bm{\delta}}^{(L,T)}({\bm{x}})-\tilde{{\bm{\delta}}}^{(L,T)}({\bm{x}})\|_{2}$
$\displaystyle\leq\sigma_{v}\|\phi^{\prime}({\bm{g}}^{(L,T)})\|_{\infty}\|({\bm{v}}-\tilde{{\bm{v}}})\|_{2}+\sigma_{v}\|{\bm{v}}\|_{2}\|\phi^{\prime}({\bm{g}}^{(L,T)}({\bm{x}}))-\phi^{\prime}({\bm{g}}^{(L,T)}({\bm{x}}))\|_{2}$
$\displaystyle\leq(\sigma_{v}C+2\sigma_{v}Ck)\|\theta-\tilde{\theta}\|_{2}\sqrt{n}.$

And similarly using same techniques we have

$\displaystyle\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})-\tilde{{\bm{\delta}}}^{(\ell,t)}({\bm{x}})\|_{2}\leq(\sigma_{w}C+3\sigma_{w}Ck+\sigma_{u}C+3\sigma_{u}Ck)\|\theta-\tilde{\theta}\|_{2}\sqrt{n}.$

As a result, there exists $K_{1}$ that is a function of
$\sigma_{w},\sigma_{u},\sigma_{b},L,T$ and the norm of the inputs.

Now we prove the local Lipchitzness of the Jacobian

$\displaystyle\|{\bm{J}}(\theta,{\bm{x}})\|_{F}\leq$
$\displaystyle\sum_{\ell=2}^{L}\sum_{t=1}^{T}\bigg{(}\frac{1}{n}\left\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})\left(\sigma_{w}^{\ell}{\bm{h}}^{(\ell,t-1)}({\bm{x}})\right)^{\top}\right\|_{F}$
$\displaystyle+\frac{1}{n}\left\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})\left(\sigma_{u}^{\ell}{\bm{h}}^{(\ell,t-1)}({\bm{x}})\right)^{\top}\right\|_{F}+\frac{1}{\sqrt{n}}\left\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})\bm{\cdot}\sigma_{b}^{\ell}\right\|_{F}\bigg{)}$
$\displaystyle+\sum_{t=1}^{T}\bigg{(}\frac{1}{n}\left\|{\bm{\delta}}^{(1,t)}({\bm{x}})\left(\sigma_{w}^{1}{\bm{h}}^{(1,t-1)}({\bm{x}})\right)^{\top}\right\|_{F}$
$\displaystyle+\frac{1}{\sqrt{nm}}\left\|{\bm{\delta}}^{(1,t-1)}({\bm{x}})\left(\sigma_{u}^{1}{\bm{x}}_{t}\right)^{\top}\right\|_{F}+\frac{1}{\sqrt{n}}\left\|{\bm{\delta}}^{(1,t)}({\bm{x}})\bm{\cdot}\sigma_{b}^{1}\right\|_{F}\bigg{)}+\frac{\sigma_{v}}{\sqrt{n}}\|{\bm{h}}^{(L,T)}({\bm{x}})\|_{F}$
$\displaystyle\leq\bigg{(}\sum_{\ell=2}^{L}\sum_{t=1}^{T}(K_{1}^{2}C\sigma_{w}^{\ell}+K_{1}^{2}C\sigma_{u}^{\ell}+\sigma_{b}^{\ell}K_{1})$
$\displaystyle+\sum_{t=1}^{T}(K_{1}^{2}C\sigma_{w}^{1}+\frac{K_{1}\sigma_{u}^{1}}{\sqrt{m}}\|{\bm{x}}_{t}\|_{2}+\sigma_{b}^{1}K_{1})+\sigma_{v}CK_{1}\bigg{)}.$

And for $\theta,\tilde{\theta}\in B(\theta_{0},R)$ we have

$\displaystyle\|{\bm{J}}(\theta,{\bm{x}})-\tilde{{\bm{J}}}(\theta,{\bm{x}})\|_{F}\leq$
$\displaystyle\sum_{\ell=2}^{L}\sum_{t=1}^{T}\bigg{(}\frac{1}{n}\left\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})\left(\sigma_{w}^{\ell}{\bm{h}}^{(\ell,t-1)}({\bm{x}})\right)^{\top}-\tilde{{\bm{\delta}}}^{(\ell,t)}({\bm{x}})\left(\sigma_{w}^{\ell}\tilde{{\bm{h}}}^{(\ell,t-1)}({\bm{x}})\right)^{\top}\right\|_{F}$
$\displaystyle+\frac{1}{n}\left\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})\left(\sigma_{u}^{\ell}{\bm{h}}^{(\ell,t-1)}({\bm{x}})\right)^{\top}-\tilde{{\bm{\delta}}}^{(\ell,t)}({\bm{x}})\left(\sigma_{u}^{\ell}\tilde{{\bm{h}}}^{(\ell,t-1)}({\bm{x}})\right)^{\top}\right\|_{F}$
$\displaystyle+\frac{1}{\sqrt{n}}\left\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})\bm{\cdot}\sigma_{b}^{\ell}-\tilde{{\bm{\delta}}}^{(\ell,t)}({\bm{x}})\bm{\cdot}\sigma_{b}^{\ell}\right\|_{F}$
$\displaystyle+\sum_{t=1}^{T}\bigg{(}\frac{1}{n}\left\|{\bm{\delta}}^{(1,t)}({\bm{x}})\left(\sigma_{w}^{1}{\bm{h}}^{(1,t-1)}({\bm{x}})\right)^{\top}-\tilde{{\bm{\delta}}}^{(1,t)}({\bm{x}})\left(\sigma_{w}^{1}\tilde{{\bm{h}}^{(1,t-1)}}({\bm{x}})\right)^{\top}\right\|_{F}$
$\displaystyle+\frac{1}{\sqrt{nm}}\left\|{\bm{\delta}}^{(1,t-1)}({\bm{x}})\left(\sigma_{u}^{1}{\bm{x}}_{t}\right)^{\top}-\tilde{{\bm{\delta}}}^{(1,t-1)}({\bm{x}})\left(\sigma_{u}^{1}{\bm{x}}_{t}\right)^{\top}\right\|_{F}$
$\displaystyle+\frac{1}{\sqrt{n}}\left\|{\bm{\delta}}^{(\ell,t)}({\bm{x}})\bm{\cdot}\sigma_{b}^{\ell}-\tilde{{\bm{\delta}}}^{(\ell,t)}({\bm{x}})\bm{\cdot}\sigma_{b}^{\ell}\right\|\bigg{)}+\frac{\sigma_{v}}{\sqrt{n}}\|{\bm{h}}^{(L,T)}({\bm{x}})-\tilde{{\bm{h}}^{(L,T)}({\bm{x}})}\|_{F}$
$\displaystyle\leq\bigg{(}\sum_{\ell=2}^{L}\sum_{t=1}^{T}(4K_{1}^{2}C\sigma_{w}^{\ell}+4K_{1}^{2}C\sigma_{u}^{\ell}+\sigma_{b}^{\ell}K_{1})$
$\displaystyle+\sum_{t=1}^{T}(4K_{1}^{2}C\sigma_{w}^{1}+\frac{K_{1}\sigma_{u}^{1}}{\sqrt{m}}\|{\bm{x}}_{t}\|_{2}+\sigma_{b}^{1}K_{1})+\sigma_{v}CK_{1}\bigg{)}\|\theta-\tilde{\theta}\|_{2}.$

The above proof can be generalized to the entire dataset by a straightforward
application of the union bound. This concludes the proof for Theorem 2.