Coupled Oscillatory Recurrent Neural Network (coRNN): An accurate and (gradient) stable architecture for learning long time dependencies
Authors: T. Konstantin Rusch
Siddhartha Mishra
Venue: ICLR 2021 Oral (2021)
Keywords: RNNs
Oscillators
Gradient stability
Long-term dependencies

Abstract: Circuits of biological neurons, such as in the functional parts of the brain can be modeled as networks of coupled oscillators. Inspired by the ability of these systems to express a rich set of outputs while keeping (gradients of) state variables bounded, we propose a novel architecture for recurrent neural networks. Our proposed RNN is based on a time-discretization of a system of second-order ordinary differential equations, modeling networks of controlled nonlinear oscillators. We prove precise bounds on the gradients of the hidden states, leading to the mitigation of the exploding and vanishing gradient problem for this RNN. Experiments show that the proposed RNN is comparable in performance to the state of the art on a variety of benchmarks, demonstrating the potential of this architecture to provide stable and accurate RNNs for processing complex sequential data.

Published as a conference paper at ICLR 2021
COUPLED OSCILLATORY RECURRENT NEURAL NET- WORK (CORNN): AN ACCURATE AND (GRADIENT)
STABLE ARCHITECTURE FOR LEARNING LONG TIME DEPENDENCIES
T. Konstantin Rusch
Siddhartha Mishra
Seminar for Applied Mathematics (SAM) Department of Mathematics
Seminar for Applied Mathematics (SAM)
Department of Mathematics
ETH Zürich Zürich, 8092, Switzerland
ETH Zürich
Zürich, 8092, Switzerland smi shraßethz. ch
t trusch@ethz.co
ABSTRACT
Circuits of biological neurons, such as in the functional parts of the brain can be modeled as networks of coupled oscillators. Inspired by the ability of these systems express a rich set of outputs while keeping (gradients of) state variables bounded, we propose a novel architecture for recurrent neural networks. Our proposed RNN is based fime-discretization of system of second-order ordinary differential equations, modeling networks of controlled nonlinear oscillators. We prove precise bounds on the gradients of the hidden states, leading the mitigation of the exploding and vamishing gradient problem for this RNN. Experiments show that the proposed RNN comparable performance the state of the art on variety of benchmarks, demonstrating the potential of this architecture to provide stable and accurate RNNs for processing complex sequential data.
INTRODUCTION
Recurrent neural networks (RNNs) have achieved tremendous success in variety of tasks involving sequential (time series) inputs and outputs, ranging from speech recognition computer vision and natural language processing, among others. However, it is well known that training RNNs to process inputs over long time scales (input sequences) is notoriously hard on account of the so-called exploding and vanishing gradient problem (EVGP) (Pascanu et al. 2013), which stems from the fact that the sell-establishe BPTT algorithm for training RNNs requires computing products of gradients (Jacobians) of the underlying hidden states very long time scales. Consequently, the overall gradient can grow (to infinity) or decay (to zero) exponentially fast with respect to the number of recurrent interactions.
variety of approaches have been suggested to mitigate the exploding and vanishing gradient problem. These include adding gating mechanisms the RNN in order to control the flow of information in the network, leading to architectures such as long short-term memory (LSTM) (Hochreiter & Schmidhuber 1997) and gated recurring units (GRU) (Cho et al., 2014). that can overcome the vanishing gradient problem on account of the underlying additive structure. However, the gradients might still explode and learning very long term dependencies remains a challenge (Li al. 2018). Another popular approach for handling the EVGP is to constrain the structure of underlying recurrent weight matrices by requiring them to be orthogonal (unitary), leading to the so-called orthogonal RNNs (Henaff et al. 2016 [ Arjovsky al. 2016) Wisdom et al. 2016 Kerg et al. (2019) and references therein. By construction, the resulting Jacobians have eigen- and singular-spectra with unit norm, alleviating the EVGP. However as pointed out by Kerg al. (201). imposing such constraints the recurrent matrices may lead to significant loss of expressivity of the RNN resulting in inadequate performance on realistic tasks.
In this article, we adopt a different approach, based on observation that coupled networks of controlled non-linear forced and damped oscillators, that arise in many physical, engineering and biological

---

Published as a conference paper at ICLR 2021
systems, such as networks of biological neurons, do seem to ensure expressive representations while constraining the dynamics of state variables and their gradients. This motivates us to propose a novel architecture for RNNs, based on time-disse deeeeretizzationss of second-order systems of non-linear ordinary differential equations (ODEs) (T) that model coupled oscillators. Under verifiable hypotheses, we able to rigorously prove precise bounds on the hidden states of these RNNs and their gradients, enabling a possible solution of the exploding and vanishing gradient problem, while demonstrating through benchmark numerical experiments, that the resulting system still retains sufficient expressivity, i.e ability to process complex inputs, with a competitive performance, with respect to the state of the art, on a variety of sequential learning tasks.
THE PROPOSED RNN
Our proposed RNN is based on the following second-order system of ODEs,
y" = a (Wy + Wy' + Vu + b) ny -cy.
(1)
Here, t E [0, 1] is the (continuous) time variable, u = u(t) E Rd is the time-dependent input signal,
y = yt) E I is the hidden state of the RNN with W, W E mm V E Rm are weight matrices, b E R is the bias vector and 0 < 7, are parameters, representing oscillation frequency and the amount of damping (friction) in the system, respectively. R R is the activation function, set to (u) = tanh(u) here. By introducing the so-called velocity variable z = y'(t) E R""
we rewrite ①1 the first-order system:
y' = 2, z' = (Wy + Wz + Vu b)
- -y
(2)
We fix a timestep 0 < At < 1 and define our proposed RNN hidden states at time t, = nAt E [0. 1] (while omitting the affine output state) as the following IMEX (implicit-explicit) discretization of the first order system 2
Yn Yn-1 + Atz,:
(3)
2, = - + Ato (Wyn-1 + Wzn + Vu, + b) Atryn-1 - Ateza:
with either π n or ń in the implicit (ñ = n) or explicit (n = n 1) treatment of the damping term -c2 in whereas both versions retain the implicit treatment of the first equation in (2).
1. Note that the only difference in the two versions of the RNN 3 lies
Motivation and background. To see that the underlying ODE (2) models a coupled network of controlled forced and damped nonlinear oscillators, we start with the single neuron (scalar) case by setting d = m = 1 in (1) and assume an identity activation function (z) x. Setting V = b- = 0 leads the simple ODE, y" + -yy== = O, which exactly models simple
W
W W- harmonic motion with frequency y, for instance that of a mass attached to a spring (Guckenheimer & Holmes| 1990). Letting e 0 in (T) adds damping or friction to the system (Guckenheimer & Holmes] 1990). Then, by introducing non-zero V in I we drive the system with driving force proportional to the input signal u(t). The parameters V modulate the effect of the driving force, W controls the frequency of oscillations and W the amount of damping in the system. Finally, the tanh activation mediates a non-linear response in the oscillator. In the coupled network (2) with m 1, each neuron updates its hidden state based on the input signal well information from other neurons. The diagonal entries of W (and the scalar hyperparameter 7) control the frequency whereas the diagonal entries of W (and the hyperparameter e) determine the amount of damping for each neuron, respectively, whereas the non- non- -oidonaa entries of these matrices modulate interactions between neurons. Hence, given this behavior of the underlying ODE (2). we term the RNN (3) as a coupled oxcillatory Recurrent Neural Network (coRNN).
The dynamics of the ODE ( (and the RNN (3) for a single neuron are relatively straightforward. As we illustrate in Fig. 6 of supplementary material SM&C input signals drive the generation of (superpositions of) oscillatory fave-forms, whose amplitude and (multiple) frequencies are controlled by the tunable parameters W, W. V, b. Adding tanh activation does not change these dynamics much. This is in contrast to truncating tanh to leading non-linear order by setting o = x x² /3. which yields Duffing type oscillator that is characterized by chaotic behavior (Guckenheimer & Holmes] 1990). Adding interactions between neurons leads to further accentuation of this generation of superposed wave forms (see Fig. 6 in SM&C) and even with very simple network topologies, one
2

---

Published as a conference paper at ICLR 2021
sees the emergence of non-trivial non-oscillatory hidden states from oscillatory inputs. In practice, a network of a large number of neurons is used and can lead to extremely rich global dynamics. Hence, we argue that the ability of network of (forced, driven) oscillators to access a very rich set of output states may lead to high expressivity of the system, allowing it to approximate outputs from complicated sequential inputs.
Oscillator networks are ubiquitous in nature and in engineering systems (Guckenheimer & Holmes 1990, Strogatz 2015) with canonical examples being pendulums (classical mechanics), business cycles (economics), heartbeat (biology) for single oscillators and electrical circuits for networks of oscillators. Our motivating examples arise in neurobiology. where individual biological neurons can be viewed as oscillators with periodic spiking and firing of the action potential. Moreover, functional circuits of the brain, such as cortical columns and increasingly interpreted by networks of oscillatory neurons, see [Stiefel & Ermentrout (2016) for an overview. Following well-established paths in machine learning, such as for convolutional neural networks (LeCun et al., 2015). our focus here is to abstract the essence of functional brain circuits being networks of oscillators and design an RNN based on much simpler mechanistic systems, such as those modeled by 2), while ignoring the complicated biological details of neural function.
circuits, are being
Related work. There is an increasing trend of basing RNN architectures on ODEs and dynamical systems. These approaches can roughly be classified into two branches, namely RNNs based on discretized ODEs and continuous-ting
RNNs. Examples of continuous-time approaches include neural ODEs (Chen et al. 20118) with ODE-RNNS (Rubanova et al.) (019) as its recurrent extension as well as (0177) and references therein, to name just a few. We focus, however, in this article on an E-Einiiired discrete-time RNN, the proposed coRNN is derived from a discretization of the ODE I A good example for a discrete-ting time ODE-based RNNs is the so-called anti-symmetric RNN of Chang al. ) where the RNN architecture is based on stable ODE resulting from a skew-symmetric hidden weight matrix, thus constraining the stable (gradient) dynamics of the network. This approach has much in common with previously mentioned unitary/orthoggna//noo normal RNNs in constraining the structure of the adding such strong constraints might reduce expressivity of the resulting RNN and might lead to inadequate performance on complex tasks. In contrast to these approaches, our proposed coRNN does not explicitly constrain the weight matrices but relies on the dynamics of the underlying ODE (and the IMEX discretization (3), to provide gradient stability. Moreover, no gating mechanisms as in LSTMs/GRUs are used in the current version of CORNN. There is also an increasing interest in designing hybrid methods, which use a discretization of an ODE (in particular a Hamiltonian system) in order to learn the continuous representation of the data, see for instance Greyydauus et al. (2019); Chen et al. (020). Overall, our approach here differs from these papers in our use of networks of oscillators to build the RNN.
layer weight matrices. However,
RIGOROUS ANALYSIS OF THE PROPOSED RNN
An attractive feature of the underlying ODE system lies in the fact that the resulting hidden states (and their gradients) bounded (see SM for precise statements and proofs). Hence, one can expect that a suitable discretization of the ODE 2b that preserves these bounds will not have exploding gradients. We claim that one such structure preserving discretization is given by the IMEX discretization that results in the RNN (3) and proceed to derive bounds on this RNN below.
Following standard practice we set y (0) z(0) 0 and purely for the simplicity of exposition, we set the control parameters, c y = 1 and π. = n in leading to,
Yn Yn-1 + Atzn; Zn T+Ar
1 At
At THA7 ( Vn-1 An-1 Wyn-1 Wz, + Vu, + b.
(4)
Analogous results and proofs for the case where ń = n 1 and for general values of e, y are provided in F
Bounds on the hidden states. As with the underlying ODE (2. the hidden states of the RNN are bounded, i.e.
3

---

Published as a conference paper at ICLR 2021
Proposition 3.1 Let Yn, Zn be the hidden states of the RNN 44 for 1 I S N, then the hidden states satisfy the following (energy) bounds:
yn 2 S nmmt mty m.
(5)
The proof of the energy bound (5) is provided SM and a straightforward variant the proof (see E.2 yields estimate the sensitivity of the hidden states changing inputs. As with the underlying ODE (see SM this bound rules out chaotic behavior of hidden states.
Bounds hidden state gradients.
We train the RNN 3 to minimize the loss function,
Σ En, En = ly, yn112. N
(6)
with y being the underlying ground truth (training data). During training, we compute gradients of the loss function (6) with respect to the weights and biases + [W, W, V, V.bl.
20
N
gE,
AOEO.
00
N
Let Yn, be the hidden states generated by the RNN (4) We assume that the time step At << can be chosen such that,
Proposition 3.2
max (11( + (
A/|
= 2 sr.l.
(8)
1+
1+ At
Denoting ô = T4 the gradient of the loss function bounded as,
& with respect to any parameter 0 + is
30
3
(m + Y /
00
2
with Y max ISn5N
Vallo be bound the underlying training data.
Sketch of the proof. Denoting X,, = -yn Zn], we can apply the chain rule repeatedly (for instance as Pascanu et al. (013) to obtain,
0E
0 HE, OX,
Σ
(10)
00
Ox,
OX 00
ISASn
De (4)
59
Here, the notation X+Xx
36 refers to taking the partial derivative of X₁ with respect to the parameter
while keeping the other arguments constant. This quantity can be readily calculated from the structure
of the RNN (4) and is presented in the detailed proof provided SM E.3 From we can directly
compute that Dr xx yy, yn,
Repeated application of the chain rule and direct calculation with 4 yields,
ox,
II keisn OX
OX;
ax,
I AIB:-
AIC,-
(11)
0x
ax
B
C.-1
where I is the identity matrix and
B jAt (diag(o((A, 1), 8 At diag(o' (Aii))))))))))) (12)
is straightforward to calculate using the assumption 8 that B < n and |C I n 7 Using the definitions of matrix norms and obtain:
0x,
< max (1 + At B + CIC,,)l)e),, |B. + | Cc-1ll)
0X
(13)
max (1 + A + 2 2m) < +3t".
4

---

Published as a conference paper at ICLR 2021
Therefore, using
have
axm 0X
ax
I {1+ (1+35))** 1+ -)
(14)
0X;
k<išn
Note that we have used an expansion around 1 and neglected terms of ( as At << 1. We remark that the bound (33 is the crux of our argument about gradient control as structure of the RNN that the recurrent matrices have close to unit norm. The detailed proof is
from the
presented in SMSE.3 As the entire gradient of the loss function (6). with respect the weights and
biases of the network, is bounded above in (9) the exploding gradient problem is mitigated for this RNN.
On the vanishing gradient problem. The vanishing gradient problem (Pascanu et al. 2013) arises if
(E)
defined in (00). - 0 exponentially fast in k, for l: << (long-term dependencies). In that case, the RNN does not have long-term memory, as the contribution of the k-th hidden state to error at time step is infinitesimally small. We already from 14 that
20
8Xx 1oo (independently of k). Thus, we should not expect the products in (10) to decay fast. In fact, we will provide much more precise characterization of this gradient. To this end, we introduce the following onder-notation,
0x
8 (a), for a, R if there exists constants C,C such that Co S Ca.
(15)
M = (a),, for M E Redi a E RR+, there exists constant C such that |M S Co.
For simplicity of notation, will also set yn U, = O. for all b 0 and
1 in
8) and we
will only consider 8 = Ws, for some 1 I j.j S m in the following proposition.
Proposition 3.3 Let n. be the hidden states generated by the RNN 4 Under the assumption that y² = O(Vt,) for all < S m and (8), the gradient for long-term dependencies satisfies,
ac
00 o ( + (25(1 (0(1+8)JA0 (15(1+ +))(a) + O(Ai)), è sech? (VKDI(I+AA)
+ +di)). k << n.
(16)
This precise bound (16) on the gradient shows that although the gradient can small, i.c O(At)) 8
is in fact independent of k. ensuring that long-term dependencies contribute to gradients at much later steps and mitigating the vanishing gradient problem. The detailed proof is presented in SM
that the RNN 33 indeed satisfied similar bounds to the underlying ODE that resulted in upper bounds on the hidden states and its gradients. However, the lower bound on the gradient (16) due the specific choice of this discretization and does not appear have continuous analogue, making the specific choice of discretization of (2) crucial for mitigating the vanishing gradient problem.
Summarizing,
EXPERIMENTS
as this version We present results on variety of learning tasks with coRNN 5 with Details of the training resulted in marginally better performance than the version with π procedure for each experiment can be found in SMęB We wish to clarify here that we straightforward hyperparameter tuning protocol based on a validation set and do not use additional performance enhancing tools, such as dropout (Srivastava al. 2014). gradient clipping (Pascanu et al. (013) or batch normalization (loffe & Szegedy 2015), which might further improve the performance of coRNNs.
Adding problem. We start with the well-known adding problem (Hochreiter & Schmidhuber 1997). proposed to test the ability of an RNN to learn (very) long-term dependencies. The input is wwo-dimensional sequence of length T. with the first dimension consisting of random numbers drawn from ((0, 11) and with two non non-zero zero entries (both set to 1) in the second dimension, chosen at random locations, but one each in both halves of the sequence. The output is the sum of two numbers
5

---

Published as a conference paper at ICLR 2021
613
612
0.122
8.100
FaSIRNN
and sym. RNN
Tarhh RNN
Eaunlite
Traiiig (tundrado)
(hundreds)
(handreds)
Figure 1: Results of the adding problem for coRNN, expRNN, FastRNN, anti.sym. RNN and tanh RNN based on three different sequence lengths T, i.c. T 500, T 2000 and 5000.
of the first dimension at positions, corresponding to the two entries in the second dimension. We compare the proposed coRNN to three recently proposed RNNs, which explicitly designed to learn LTDs, namely the FastRNN (Kusupati et al. (018), the antisymmetric (anti.sym.) RNN (Chang) & Martínez-Rubil lanh RNN, with the goal of beating the baseline mean square error (MSE) of 0.167 (which stems from the variance of the baseline output 1). All methods have 128 hidden units (dimensionality of the hidden state y) and the same training protocol is used in all cases. Fig. T shows the results for different lengths 7 of the input sequences. We can see that while the tanh RNN not able beat the baseline for any sequence length, the other methods successfully learn the adding task for T = 500. However, in this case, coRNN converges significantly faster and reaches a lower test MSE than other tested methods. When setting the length the much more challenging case of T 2000, that only CoRNN and the expRNN beat the baseline. However, the expRNN fails to reach desired test MSE of 0.01 within training time. In order further demonstrate the superiority of coRNN over recently proposed RNN architectures for learning LTDs, we consider the adding problem for T 5000 and observe that coRNN converges very quickly even in this case, while expRNN fails consistently beat the baseline. We thus conclude that the coRNN mitigates the vanishing/exploding gradient problem even for very long sequences.
al. 2019) and the expRNN (Lezcano-Casate
2019), and plain vanilla
Table I: Test accuracies on sMNIST and PSMNIST (we provide our own psMNIST result for the FastGRNN, official result for this task has been published far).
Model
sMNIST
PSMNIST
units
# params
uRNN
Arjovsky et al. 2016
95.1%
91.4%
512
9k
LSTM Helfrich et 2018
98.9%
92.9%
256
270k
GRU Chang et al. 2017)
99.1%
94.1%
256
200k
anti.sym. RNN (Chang et al.
2019
98.0%
95.8%
128
10k
DTRIVoo (Casado 2019
99.0%
96.8%
512
137k
FastGRNN (Kusupati et al. 2018)
98.7%
94.8%
128
18k
coRNN (128 units)
99.3%
96.6%
128
34k
coRNN (256 units)
99.4% %
97.3%
256
134k
2015) is a benchmark Sequential (permuted) MNIST. Sequential MNIST (SMNIST) (Le et al.
for RNNs, in which the model is required to classify an MNIST (LeCun et al. 1998 digit one pixel at time leading to a classification task with a sequence length of 784. In permuted sequential MNIST (psMNIST), a fixed random permutation is applied in order to increase the time-delay between interdependent pixels and to make the problem harder. In Table []] we compare the test accuracy for CoRNN on sMNIST and PsMNIST with recently published best case results for other recurrent models, which were explicitly designed to solve long-term dependencies together with baselines corresponding to gated and unitary RNNs. To the best of our knowledge the proposed coRNN outperforms all single-layer recurrent architectures, published in the literature, for both the sMNIST and PsMNIST. Moreover in Fig.2) we present the performance (with respect to number of epochs) of different RNN architectures for PSMNIST with the same fixed random permutation and the

---

Published as a conference paper at ICLR 2021
number of hidden units, ie. 128. As seen from this figure, coRNN clearly outperforms the other architectures, some which were explicitly designed to learn LTDs, handily for this permutation.
1.0
0.10
505
B Adt
line:
585
cORMN
DOMNES
raaay
dashed:
580
575
3007
FarGGRNN
0/44
#MNIST
570
0.6
50.0
5000
35000
epoch
Learing steps
Figure Performance on PsM- NIST for different models, all with 128 hidden units and the same fixed random permutation.
Figure 3: Weight assumptions
Figure 4: Ablation study on the hyperparameters c in 33
evaluated dur- ing training for all LTD experi- ments (mean and standard devia- experiment. tion of 10 different runs for each task).
(8), with
ing the noise padded CIFAR-10
Noise padded CIFAR-I0. Another challenging test problem for learning LTDs is the recently proposed noise padded CIFAR-10 experiment by [Chang et al. (2019), in which CIFAR- 10 data points (Krizhevsky et al. 2009) are fed the RNN row-wise and flattened along the channels resulting in sequences of length 32. To test the long term memory, entries of uniform random numbers are added such that the resulting sequences have length of 1000, i.e. the last 968 entries of each sequence are only noise to distract the network. Table 2 shows the result for coRNN together with other recently published best case results. We observe that CoRNN readily outperforms other RNN architectures on this benchmark, while requiring only 128 hidden units.
Table 2: Test accuracies noise padded CIFAR-10.
Model
test accuracy
# units
params
LSTM (Kag et al.
2020
11.6%
128
64k
Incremental RNN Kag al. 2020
54.5%
128
12k
FastRNN (Kag et al. 2020)
45.8%
128
16k
anti.sym. RNN (Chang
2019
48.3%
256
36k
Gated anti syti.sym. RNN (Chang et al.
2019
54.7%
256
37k
Lipschitz RNN Erichson 2020
55.2%
256
134k
coRNN
59.0%
128
46k
Human activity recognition. This experiment is based on the human activity recognition data set provided by |Aggutaa al.) (2012). The data set is collection of tracked human activities, which were measured by an accelerometer and gyroscope on a Samsung Galaxy S3 smartphone. Six activities were binarized to obtain two merged classes (Sitting, Laying, Walking Valking__psttirrs) and (Standing, Walking, Walking_Downsteeee Downstairs), leading to the HAR-2 data set, which was first proposed in Kusupati et al. (2018). Tshos the result for coRNN together with other very recently published best case results on the same data set. We can see that CoRNN readily outperforms all other methods. We also ran this experiment even in this case, the tiny coRNN beats all baselines. We thus conclude that coRNN can efficiently be used on constrieed IoT cono-coontsllers.
tiny coRNN with very few parameters, i.e. only 1k. We can that
IMDB sentiment analysis. The IMDB data set (Maas et al. 2011 is collection of 50k movie reviews, where 25k reviews are used for training (with 7.5k of these reviews used for validating) and 25k reviews used for testing. The aim of this binary sentiment classification task is decide whether a movie review is positive or negative. We follow the standard procedure by initializing the word embedding with pretrained 100d GloVe (Pennington al. 2014) vectors and restrict the
7

---

Published as a conference paper at ICLR 2021
Table Test accuracies on HAR-2.
Model
test accuracy
# units
# params
GRU (Kssupati et al. 2018
93.6%
75
19k
LSTM (Kag et al. 2020)
93.7%
64
16k
FastRNN Kusupati et al.
2018
94.5%
80
7k
FastGRNN (Kusupati a.. 2018
95.6%
80
anti.sym. RNN (Kag ct al. 2020)
93.2%
120
8k
incremental RNN
(Kag ct al.
2020
96.3%
64
4k
coRNN
97.2%
64
9k
tiny CORNN
96.5%
20
1k
dictionary 25k words. Table/ssowss shows the results for coRNN and other recently published models, which are trained similarly and have the same number of hidden units, i.e. 128. We can see that coRNN compares favorable with gated baselines (which are known perform very well on this task), while at the same time requiring significantly less parameters.
Table 4: Test accuracies on IMDB.
Model
test accuracy units
params
LSTM (Campos et al. 2018)
86.8%
128
220k
Skip LSTM(Campos al.
2018
86.6%
128
220k
GRU (Campos et al. 20188
86.2%
128
164k
Skip GRU (Campos al. 2018
86.6%
128
164k
ReLU GRU (Dey & Salemt 2017)
84.8%
128
99k
coRNN
87.4%
128
46k
Further experimental results. To shed further light the performance of coRNN, consider the following issues. First, the theory suggested that coRNN mitigates the exploding/vanishing gradient problem as long as the assumptions (8) on the time step At and weight matrices W, W hold.
Clearly one can choose suitable At to enforce 8 before training, but do these assumptions remain valid during training? In SM$E.4 we argue, based on worst-case estimates, that the assumptions will remain valid for possibly large number of training steps. More pertinently, we verify experimentally that (8) holds during training. This demonstrated Fig. 3 where we show that (8) holds for all LTD tasks during training. Thus, the presented theory applies and one can expect control over hidden state gradients with coRNN. Next, we recall that the frequency parameter 7 and damping parameter play role for coRNNs (see for the theoretical dependence and Table 8 for best performing values of y for each numerical experiment within the range considered in Table (7). How sensitive is the performance of coRNN to the choice of these 2 parameters? To investigate this dependence, we focus on the noise padded CIFAR-10 experiment and show the results of an ablation study in Fig. 4 where the test accuracy for different CoRNNs based on
dimensional hyperparameter grid (e, 7) E [0.8, 1.8) x [5.7, 17, 7] (i.e., sufficiently large intervals around the best performing values of e, from Table 8) is plotted. We observe from the figure that although there are reductions in test accuracy for non- optimall values of (c, Y), there is no large variation and the performance is rather robust with respect these hyperparameters. Finally, note that we follow standard practice and present best reported results with coRNN well as other competing RNNs in order to compare the relative performance. However, natural investigate the dependence of these best results on the random initial (before training) values of the weight matrices. To this end, in Table SM, we report the mean and standard deviation (over 10 retrainings) of the test accuracy with coRNN various learning tasks and find that the mean value is comparable the best reported value, with low standard deviations. This indicates further robustness of the performance of coRNNs.
8

---

Published as a conference paper at ICLR 2021
Table 5: Distributional information (mean and standard deviation) on the results for each classification experiment presented in the paper based on 10 re-trainings of the best performing coRNN using random initialization of the trainable parameters.
Experiment
Mean
Standard deviation
sMNIST (256 units)
99.17%
0.07%
PSMNIST (256 units)
96.10%
1.20%
Noise padded CIFAR-10
58.56%
0.35%
HAR-2 (64 units)
96.01%
0.53%
IMDB
86.65%
0.31% %
DISCUSSION
Inspired by many models in physics, biology and engineering, we proposed a novel RNN architecture (3) based on a model ①T of network of controlled forced and damped oscillators. For this RNN, we rigorously showed that under verifiable hypotheses on the time step and weight matrices, the hidden states are bounded (5) and obtained precise bounds on the gradients (Jacobians) of the hidden states, (9) and (1)). Thus by design, this architecture can mitigate the exploding and vanishing gradient problem (EVGP) for RNNs. We present series of numerical experiments that include sequential image classification, activity recognition and sentiment analysis, to demonstrate that the proposed coRNN keeps hidden states and their gradients under control, while retaining sufficient expressivity to perform complex tasks. Thus, we provide a novel and promising strategy for designing RNN architectures that motivated by the functioning of natural systems, have rigorous bounds hidden state gradients and are robust, accurate, straightforward to train and cheap to evaluate.
built with This work can be extended in different directions. For instance in this article, we have mainly focused on the learning of tasks with long-term dependencies and observed that coRNNs are comparable in performance to the best published results in the literature. Given that coRNNs networks of oscillators, it is natural to expect that they will perform very well on tasks with oscillatory inputs/outputs, such the time series analysis of high-resolution biomedical data, for instance and EMG (electromyography data and seismic activity data from geoscience. This will be pursued in follow-up up article. Similarly, applications of coRNN to language modeling will be covered in future work.
EEG (electroencephaloggaphh
However, it essential to point out that coRNNs might not be suitable for every learning task involving sequential inputs/outputs. As concrete example, we consider the problem of predicting time series corresponding a chaotic dynamical system. We recall that by construction, the underlying ODE (2) (and the discretization (3) do not allow for super-linear (in time) separation of trajectories for nearby inputs. Thus, we cannot expect that coRNNs will be effective at predicting chaotic time series and
is indeed investigated and demonstrated for Lorenz- 96 ODE in SMŞA where observe that the coRNN is outperformed by LSTMs in the chaotic regime.
Our main theoretical focus in this paper was to demonstrate the possible mitigation of the exploding and vanishing gradient problem. On the other hand, we only provided some heuristics and numerical evidence why the proposed RNN still has sufficient expressivity. A priori, it is natural think that the proposed RNN architecture might introduce a strong bias towards oscillatory functions. However, as we argue in SM$C the proposed coRNN can be significantly more expressive, as the damping, forcing and coupling of several oscillators modulates nonlinear response to yield a very rich and diverse set of output states. This is also evidenced by the ability of coRNNs to deal with many tasks in our numerical experiments, which do not have an explicit oscillatory structure. This sets the stage for a rigorous investigation of universality of the proposed CoRNN architecture, as in the case of echo state networks in Grigoryeva & Ortega (0088). A possible approach would be to leverage the ability of the proposed RNN to convert general inputs into a rich set of superpositions of harmonics (oscillatory wave forms). Moreover, the proposed RNN was based on the simplest model ([). Much more detailed models of oscillators are available, particularly those that arise in the modeling of biological neurons, Stiefel & Ermentrout (2016) and references therein. An interesting variant of our proposed RNN would be to base the RNN architecture on these more elaborate models, resulting in analogues of the spiking neurons model of Maass (200)) for RNNs.
of coupled oscillators

---

Published as a conference paper at ICLR 2021
REFERENCES
Davide Anguita, Alessandro Ghio, Luca Oneto, Xavier Parra, and Jorge I. Reyes-Ortiz. Human activity recognition on smartphones using a multiclass hardware-friende friendly support vector machine. In International Workshop on Ambient Assisted Living, pp. 216-223. Springer, 2012.
Martin Arjovsky, Amar Shah, and Yoshua Bengio. Unitary evolution recurrent neural networks. In International Conference Machine Learning, pp. 1110-1122 1128, 2016.
Victor Campos, Brendan Jou, Xavier Giró-i-Nieto, Jordi Torres, and Shih-Fu Chang. Skip RNN: learning skip state updates in recurrent neural networks. In 6th International Conference Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 May 3, 2018, Conference Track Proceedings, 2018.
Mario Lezcano Casado. Trivializations for gradient-based optimization on manifolds. In Advances in Neural Information Processing Systems, pp. 9154-9164 2019.
Bo Chang, Minmin Chen, Eldad Haber, and Ed H. Chi. Antisymmetriemm: A dynamical system view on recurrent neural networks. 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019, 2019.
Shiyu Chang, Yang Zhang, Wei Han, Mo Yu, Xiaoxiao Guo, Wei Tan, Xiaodong Cui, Michael Witbrock, Mark A Hasegawa-Johnson, and Thomas S Huang. Dilated recurrent neural networks. In Advances in Neural Information Processing Systems, pp. 77 77-87, 20177,2017.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, pp. 6571-6583, 2018.
Zhengdao Chen, Jianyu Zhang, Martín Arjovsky, and Léon Bottou. Symplectic recurrent neural networks. In 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020, 2020.
Kyunghyun Cho, B van Merrienboer, Caglar Gulcehre, F Bougares, Schwenk, and Yoshua Bengio. Learning phrase representations using rnn encoder-decoder for statistical machine translation. In Conference on Empirical Methods in Natural Language Processing (EMNLP 2014), 2014.
Rahul Dey and Fathi M Salemt. Gate-variants of gated recurrent unit (gru) neural networks. In 2017 IEEE 60th International Midwest Symposium on Circuits and Systems (MWSCAS), pp. 1597-1600. IEEE, 2017.
Weinan E. A proposal on machine learning via dynamical systems. Commun. Math. Stat, 5:1-11, 2017.
N Benjamin Erichson, Omri Azencot, Alejandro Queiruga, and Michael W Mahoney. Lipschitz recurrent neural networks. arXiv preprini arXiv:2006. 00:::20006//007, 2020.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In Advances in Neural Information Processing Systems, pp. 15389,15389, 2019.
Lyudmila Grigoryeva and Juan-Pablo Ortega. Echo state networks are universal. Neural Networks, 108:495 508, 2018. ISSN 0893-6080.
J. Guckenheimer and P. Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer Verlag, New York, 1990.
S. Shinomoto H. Sakaguchi and Y. Kuramoto. Local and global self-entrainment in oscillator lattices. Progress of Theoretical Physics, 77:1005-1010, 1987.
Kyle Helfrich, Devin Willmott, and Qiang Ye. Orthogonal recurrent neural networks with scaled cayley transform. In International Conference on Machine Learning, pp. 1969-1978. PMLR, 2018.
10

---

Published as a conference paper at ICLR 2021
Mikael Henaff, Arthur Szlam, and Yann LeCun. Recurrent orthogonal networks and mong-memrry tasks. In Maria Florina Balcan and Kilian Q. Weinberger (, Proceedings of The 33rd Interna- tional Conference Machine Learning, volume 48 of Proceedings Machine Learning Research, pp. 2034-2042 2016.
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8):
1735-1780, 1997.
Sergey loffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Proceedings of the 32nd International Conference on Machine Learning, ICMI., volume 37 of JMLR Workshop and Conference Proceedings, pp- 448 456.
Jog,.org, 2015.
Anil Kag, Ziming Zhang, and Venkatesh Saligrama. Rnns incrementally evolving equilibrium manifold: A panacea for vanishing and exploding gradients? In 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020, 2020.
Giancarlo Kerg, Kyle Goyette, Maximilian Puelma Touzel, Gauthier Gidel, Eugene Vorontsov, Yoshua Bengio, and Guillaume Lajoie. Non-normal recurrent neural network (nnrnn); learning long time dependencies while improving expressivity with transient dynamics. In Advances in Neural Information Processing Systems, pp. 33591-13601, 2019.
Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
Aditya Kusupati, Manish Singh, Kush Bhatia, Ashish Kumar, Prateck Jain, and Manik Varma. Fastgrnn: A fast, accurate, stable and tiny kilobyte sized gated recurrent neural network. In Advances in Neural Information Processing Systems, pp. 9017-9028, 2018.
Thomas Laurent and James von Brecht. A recurrent neural network without chaos. In 5th Interna- tional Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net. 2017.
Quoc V Le, Navdeep Jaitly, and Geoffrey E. Hinton. A simple way to initialize recurrent networks of rectified linear units. arXiv preprint arXiv: 150,000100 2015.
Y. LeCun, Y. Bengio, and G. Hinton. Deep learning. Nature, 521:436 444, 2015.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278-23) -2224, 1998.
Mario Lezcano-Casado and David Murbinez-Rubio. Cheap orthogonal constraints in neural networks: A simple parametrization of the orthogonal and unitary group. volume 97 of Proceedings of Machine Learning Research, pp. 3794-3803 Long Beach, California, USA, 09-15 Jun 2019. PMLR.
Shuai Li, Wanqing Li, Chris Cook, Ce Zhu, and Yanbo Gao. Independently recurrent neural network (indran): Building a longer and deeper rnn. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 5457-5466, 2018.
Edward N Lorenz. Predictability: A problem partly solved. Proc. Seminar on Predictubility, volume 1, 1996.
Andrew L Maas, Raymond E Daly, Peter T Pham, Dan Huang, Andrew Y Ng, and Christopher Potts. Learning word vectors for sentiment analysis. In Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies, volume 1, pp. 142-150. Association for Computational Linguistics, 2011.
W. Maass. Fast sigmoidal networks via spiking neurons. Neural Computation, 9:279-304, 2001.
Razvan Pascanu, Tomas Mikoloy, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In Proceedings of the 30th International Conference on International Conference on Machine Learning, volume 28 of ICML''1, pp. I11-1310-11--318. JMLR.org, 2013.
11

---

Published as a conference paper at ICLR 2021
Jeffrey Pennington, Richard Socher, and Christopher D Manning. Glove: representation. In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), pp. 1532-1543, 2014.
Global vectors for word
Yulia Rubanova, Ricky Q. Chen, and David Duvenaud. Latent ordinary differential equations for ireeularly-sampled time series. In Advances in Neural Information Processing Systems 32, pp. 5320-5330. 2019.
Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: simple way to prevent neural networks from overfitting. The Journal of Machine Learning Research, 15(1):1929-1958 2014.
K. M. Stiefel and G. B. Ermentrout. Neurons as oscillators. Journal of Neurophysiology, 116: 1950-29600 2016.
Strogatz. Nonlinear Dynamics and Chaos. Westview, Boulder CO, 2015.
Strogatz. Exploring complex networks. Nature, 410:268-276, 2001.
A. Winfree. Biological rhythms and the behavior populations coupled oscillators. Journal of Theoretical Biology, 16:15-42. -211-4421111 1967.
Scott Wisdom, Thomas Powers, John Hershey, Jonathan Le Roux, and Les Atlas. Full-capacity unitary recurrent neural networks. In Advances in Neural Information Processing Systems, pp. 4880-4888, 2016
12

---

Published as a conference paper at ICLR 2021
Supplementary Material for:
Coupled Oscillatory Recurrent Neural Network (coRNN): An accurate and (gradient) stable architecture for learning long time dependencies
A CHAOTIC TIME-SERIES PREDICTION.
According to proposition F. CORNN does not exhibit chaotic behavior by design. While this property highly desirable for learning long-term dependencies slight perturbation of the input should not result in unbounded perturbation of the prediction), impairs the performance tasks, where the network has to learn actual chaotic dynamics. To test this numerically, we consider following version of the Lorenz 96 system: (Lorenz, 1996):
z'j (Ii+1 2,-2 (2-2));;111 Z; F.
(17)
where z E IR for all = 1...., and F is an external force controlling the level of chaos in the system. Fig 5 shows trajectory of the system (1] plotted on the 2 22-2--laee for small external
0.9
10
1.2
1.0
6
0.8
0.6
-3
0.4
0.6
0.8
1.0
1.2
Figure Exemplary (æ1, (e,x-))jjectooiies of the Lorenz 96 system for different forces F.
force of F 0.9 well trajectory for large external force of F = 8. We can that while for F = 0.9 the system does not exhibit chaotic behavior, the dynamics for F is already highly chaotic.
Our task consists of predicting the 25-th next state of a trajectory of the system 17 We provide
128 trajectories of length 2000 for each of the training, validation and test sets. The trajectories generated by numerically solving the system (17) and evaluating it at 2000 equidistantly distributed discrete time points with distance 0.01. The initial value for each trajectory chosen uniform random on [F 1 F + 1/2 around the equilibrium point (F....,F) of the system
Since LSTMs known be able produce chaotic dynamics, even in the autonomous (zero-entry) case (Laurent & von Brecht 2017). we expect them to perform significantly better than coRNN if the underlying system exhibits strong chaotic behavior. Table shows the normalized root mean square error (NRMSE) (RMSE divided by the root mean square of the target trajectory) on the test set for coRNN and LSTM. We can see that indeed for the non-chaotic case of using an external force of F = 0.9 LSTM and coRNN perform similarly. However, when the dynamics get chaotic (in this case using an external force of 8), the L.STM clearly outperforms CoRNN.
Table Test NRMSE on the Lorenz 96 system 17 for coRNN and LSTM.
Model
F= 0.9
=
units
# params
LSTM
2.0 x 10-2
6.8 x 10-2
44
9k
coRNN
2.0 10-2
9.8 10-2
64
9k
13

---

Published as a conference paper at ICLR 2021
B TRAINING DETAILS
The IMDB task was conducted on an NVIDIA GeForce GTX 1080 Ti GPU, while all other experi- ments were run Intel Xcon E3-1585Lv5 CPU. The weights and biases of coRNN are randomly initialized according to 21( 1
where Rai denotes the input dimension of each affine transformation. Instead of treating the parameters A1, and fixed hyperparameters, also treat them as trainable network parameters by constraining At to [0, 1] by using a sigmoidal activation function and e, ? by the use of ReLU for instance. However, in this case no major difference in performance is obtained. The hyperparameters are optimized with random search algorithm, where the results of the best performing coRNN (based on the validation set) are reported. The ranges of the hyperparameters for the random search algorithm are provided in Table 7 Table shows the rounded hyperparameters the best performing coRNN architecture resulting from random search algorithm for each learning task. We used 100 training epochs for sMNIST, PSMNIST and noise padded CIFAR-10 with additional 20 epochs which the learning rate was reduced by a factor of 10. Additionally, used 100 epochs for the IMDB task and 250 epochs for the HAR-2 task.
Table 7: Setting for the hyperparameter optimization of coRNN. Intervals denote ranges of the corresponding hyperparameter for the grid search algorithm, while fixed numbers mean that no hyperparameter optimization was done in this case.
task
learning rate
batch size
At
Y
Adding
x 10-2
[10-2,10-]
1,00
[1, 100]
SMNIST (Ahid 128)
[10-4 4,10-'
120
110-2,10-1
[10-|,00
[10-1,10]
sMNIST (naidd = 256)
110",10-1
120
[10-2,10-1
[10-1,10]
[10-1,10]
PsMNIST (Aijid -128)
[10 1004,10-1
120
[10-,11011
[10-1,10
[10-1,10]
PSMNIST (Nhid -256)
10 10-4,10-1
120
10-2,10-1
10-1,10
[10-1,10]
Noise padded CIFAR-10
10 10-4,10-1
100
10-2,10-1
1,100]
1,,100
HAR-2
10-4.10-
[10 10-2,10-1
10-1,10
11011,00 10
IMDB
10 10-4,10-1
64
[0 10-2,10-1
10-1,10
10-1,10
Table 8: Rounded hyperparameters of the best performing coRNN architecture.
task
learning rate
batch size
A4
Adding (T 5000)
2 x 10-2
50
1.6 x 10-2
94.5
9.5
SMNIST (Nhid = 128)
3.5 x 10-3 3
120
5.3 x 10-2
1.7
4
SMNIST (Nhid 256)
2.1 x 10
120
4.2 10
2.7
4.7
PSMNIST (Nnid = 128)
3.7 x 10-3
120
8.3 x 10-2
1.3 x 10-1
4.1
PSMNIST (nyid = -256)
5.4 10-3
120
7.6 x 10-2
4 10-1
8.0
Noise padded CIFAR-10
7.5 x 10
100
3.4 10
1.3
12.7
HAR-2
1.7 10
64
10-1
2 10-1
6.4
IMDB
6.0 X 10-4
64
5.4 x 10
4.9
4.8
C HEURISTICS OF NETWORK FUNCTION
At the level of a single neuron, the dynamics of the RNN is relatively straightforward. We start with the scalar case, i.e. m
1 and illustrate different hidden states y as function of time, for different input signals, in Fig. 6 In this figure, we consider two different input signals, one oscillatory signal given by u(t) cos( 4f and another is a combination of step functions. First, we plot the solution yt) of (I) with the parameters V, b, W, W, 0 and > 1. This simply corresponds to the case of a simple harmonic oscillator (SHO) and the solution described by a sine wave with the natural frequency of the oscillator. Next, introduce forcing by the input signal by setting V 1 and the activation function is the identity a (x) = x, leading to a forced damped oscillator (FDO), As seen from Fig. in the of an oscillatory signal, this leads very minor change the SHO,
14

---

Published as a conference paper at ICLR 2021
whereas for the step function, the change is only in the amplitude of the wave. Next, add damping by setting с = 0.25 and see that the resulting forced damped oscillator (FDO), merely damps the amplitude of the waves, without changing their frequency. Then, consider the case of controlled oscillator (CFDO) by setting W == -2 V = 2, b = 0.25, W = 0.75. As seen from Fig. G this leads significant change in the wave form in both cases. For the oscillatory input, the output is now superposition of many different forms, with different amplitudes and frequencies (phases) whereas for the step function input, the phase is shifted, Already, we can see that for linear controlled oscillator, the output can be very complicated with the superposition of different waves. This holds true when the activation function is set (zx) tanh(z) (which is our proposed coRNN). For both inputs. the output is modulated version of the one generated by CFDO, expressed as superposition of
waves. On the other hand, we also plot the solution with a Duffing type oscillator (DUFF) by setting the activation function as,
(18) In this case, the solution is very different from the CFDO and CoRNN solutions and is heavily damped (either in the output or its derivative). On the other hand, given the chaotic nature of the dynamical system in this case, a slight change in the parameters led to the output blowing up. Thus, a bounded nonlincarity seems essential in this context.
a(x) = x
3
Coupling neurons together further accentuates this generation of superpositions of different wave- forms, as seen even with the simplest case of a network with two neurons, shown in Fig. O (Bottom row). For this figure, we consider two neurons, i.e m = 2 and two different network topologies. For the first, we only allow the first neuron to influence the second one and not vice versa. This is enforced with the weight matrices,
2 0
0.75 0
W =
W
3 2
-1
0.75
We also set V = [2,2| b= (0.25, 0.25]T. Note that in this case (we name as ORD (for ordered connections)), the output of the first neuron should be exactly the same in the uncoupled (UC) case, whereas there is a distinct change in the output of the second neuron and we see that the first neuron has modulated sharp change in the resulting output wave form. is well illustrated by the emergence of an approximation to the step function (Bottom Right of Fig. 6). even though the input signal is oscillatory.
Next, we consider the case of fully connected (FC) neurons by setting the weight matrices as,
-2 1 3 -2 W
0.75 0.3
W
-1 0.75
The resulting outputs for the first neuron are now slightly different from the uncoupled case. On the the other hand, the approximation of step function output for the second neuron is further accentuated.
Even these simple examples illustrate the functioning of a network of controlled oscillators well. The input signal is converted into a superposition of waves with different frequencies and amplitudes, with these quantities being controlled by the weights and biases in (I.. Thus, very complicated outputs can be generated by modulating the number, frequencies and amplitudes of the waves. In practice, network of a large number of neurons is used and can lead to extremely rich global dynamics, along bistable heterogeneous behavior seen in systems of
the lines of emergence of synchronization idealized oscillators and explained by their mean field limit, see Sakaguchi & Kuramoto (1987); Winfree 1967)); [Srogaatz (2001). Thus, we argue that the ability of the network of (forced, driven) oscillators to access a very rich set of output states can lead to high expressivity of the system. The training process selects the weights that modulate frequencies, phases and amplitudes of individual neurons and their interaction to guide the system to its target output.
D BOUNDS ON THE DYNAMICS OF THE ORDINARY DIFFERENTIAL EQUATION
In this section, present bounds that show how the continuous time dynamics of the ordinary differ- ential equation (2) modeling non-linear damped and forced networks of oscillators, is constrained. We start with the following estimate the energy of the solutions of the system
15

---

Published as a conference paper at ICLR 2021
24
15-
Figure 6: Illustration of the hidden state y of coRNN (3) with a scalar input signal u (Top, Middle, Left) with one neuron with state (Top and Middle, Right) and two neurons with states y1 (Bottom left), and Y2 (Bottom right), corresponding to scalar input signal, shown in Top Left. Legend is SHO (simple harmonic oscillator), FHO (forced oscillator). FDO (forced and damped oscillator), CFDO (controlled forced and damped oscillator), DUFF (Duffing type) UC (Uncoupled), Ord (ordered coupling) and FC (fully coupled). Legend explained in the text.
Proposition D.1 Let y(t), z(t be the solutions of the ODE system at any time E [0, T] and assume that the damping parameter c 2 and the initial data for is given by,
(0) z00) 0.
Then, the solutions are bounded as,
mt
y(1) y(1) z(1)=z(0) 2(1) ml, 1I E (0, Tj.
(19)
To prove this proposition, we multiply the first equation in with yt) and the second equation in with 1 (t) to obtain,
d
y()) ((t)
z(t) x((Tx(()
z(1) x(()))(A()) =2(() Tz(1)
(20)
di
2
27
7
with
(t) = Wy( - Wz(t) + Vu(t) + b.
Using the elementary Cauchy's 's inequality repeatedly in (20) results in,
d. yt) (t)
z(t) (()))(t)) 2t))
1
c
z
dt
27
2n
7
VA
m
(as |o| < 1 and e e 5)
27
Integrating the above inequality over the time interval [0, and using the fact that the initial data (0) z0) = 0, we obtain the bounds (19).
The above proposition and estimate 19 clearly demonstrate that the dynamics of the network of coupled non non-linear oscillators (I) is bounded. The fact that the nonlinear activation function o = tanh is uniformly bounded in its arguments played erucial role in deriving the energy bound (19). A straightforward adaptation of this argument leads to the following proposition about the sensitivity of the system to inputs,
16

---

Published as a conference paper at ICLR 2021
Proposition D.2 (t), z(t) be the solutions the ODE system (2) with respect to the input signal u(t). Let ( z(t) be the solutions of the ODE system (2), but with respect to the input signal u(t). Assume that the damping parameter 2 and the initial data given by.
y()) = z = (0) = z(0) = O.
Then we have the following bound,
(y(t) - v(t)) (y(t) - y(t)) - 4mt
(z(t) z(t)) T (z(t) - z(t)) & Amt, Vt E (0,
7
(21)
Thus from the bound 21 there can be atmost linear separation (in time) with respect to the trajectories of the ODE for different input signals. Hence, chaotic behavior, which is characterized by the (super-)exponential separation trajectories is ruled out by the structure of the ODE system (2. Note that this property of the ODE system was primarily result of the uniform boundedness of the activation function a. Using a different activation function such as ReLU might enable to obtain an exponential separation of trajectories that a prerequisite for chaotic dynamical system.
D.1 GRADIENT DYNAMICS FOR THE ODE SYSTEM
Let @ denote the i, j-th entry of the Weight matrices W, W. V or the i-th entry of the bias vector b. We are interested in finding out how the gradients of the hidden state y (and the auxiliary hidden state 2) with respect to parameter 0, vary with time. Note that these gradients are precisely the objects of interest in the training of an RNN, based on a discretization of the ODE system To this end, we differentiate (2) with respect to the parameter f and denote
dy 30
Oz
yo(t)
(t), zo(t) = 00 (t),
obtain,
Y/ = 20:
(22)
zg diag(o'(A)) [ yo Wea] + Zia (A)p Tyo EZg.
As introduced before, Zini (A) E m is a matrix with all elements are zero except for the (i.j)-th entry which is set to o' d'(A (t)), i.c. the i-th entry of (A)), and we have,
y,
m, if 8 W i.j
z, m = =m, f = W
P u, 11 m = 1, A = b;-
d.
if f Vijz
We see from (2) that the ODEs governing the gradients with respect to the parameter also represent a system of oscillators but with additional coupling and forcing terms, proportional to the hidden states y, z or input signal u. As we have already proved with estimate (19) that the hidden states are always bounded and the input signal is assumed to be bounded, it is natural to expect that the gradients of the states with respect following proposition, which for simplicity of exposition, we consider the case of f Wi,j, as the other values of 0 are very similar in their behavior.
are also bounded. We make this statement explicit in the
Proposition D.3 Let 0 Was and y, be the solutions of the ODE system weights and the damping parameter satisfy,
Assume that the
|| || Se
then we have the following bounds on the gradients,
(o(t) (o(t) + (zo(t) 2o < (yo(0] ( yo(0)= + (z)(0) T20(00) cot mt?
te (0,7),
272
|| |1
= max
1 + ||
(23)
17

---

Published as a conference paper at ICLR 2021
The proof of this proposition follows exactly along the same lines the proof of proposition D.I and we skip the details, while noting the crucial role played by the energy bound (19
We remark that the bound (23) indicates that as long the initial gradients with respect to are bounded and the weights are controlled by the damping parameter, the hidden state gradients remain bounded in time.
E SUPPLEMENT TO THE RIGOROUS ANALYSIS OF CORNN
In this section, we supplement the section on the rigorous analysis of the proposed RNN 44. We start with
E. PROOF OF PROPOSITION
We multiply (y.T1,, 2 to
3 and use the elementary identities,
(a- b) ala b
b b 1
2(a ba bb) ba b)
a b) (a b):
2
2
to obtain the following,
Yn Yn + 2 nn
y y2 + nn 12n11
(yn Thn-yn-ie Yn-1) (yn Vn-1)
2
in -n-
Zn-1) (Zn
+ Alz o(An-1) Alu 2n
2
y y + 2n 1
12111
AI (1/2 + A1/2 - 1) Un Un + At T (An-1))o(A,,
y nn yn- + 1n-11nn1
mAt
0² < sandd e <<1
2
Iterating the above inequality n times leads to the energy bound,
ynyn 2 nn S vo
20 20 nmAt md..
(24)
as yo Zo O.
E.2 SENSITIVITY TO INPUTS
Next, examine how changes in the input signal affect the dynamics. We have the following proposition:
Proposition E.1 Let yn: Zn be the hidden states of the trained RNN with respect the input N [un]] 1==1 T T n= then the differences in the hidden states are bounded by,
and let Fn, Zn be the hidden states of the same RNN (4) but with respect to the input
(yn -.)) (yn -x) (z,, Zn) (z,,
Zn) S 4mt
(25)
The proof of this proposition is completely analogous the proof of proposition we subtract
y, y
Atz
2 T+at + I+
o (An-)) - At I4A J. 1- A WF + WZ + Vū, + b.
At
(26)
from and multiply (vyn (z,, Zn,7) the difference. The estimate (5) follows
identically to the proof of (5) (presented above) by realizing that (An-1) (An-1) S 2.
Note that the bound 25 ensures that the hidden states can only separate linearly in time for changes in the input. Thus, chaotic behavior, such as for Duffing type oscillators, characterized by at least exponential separation of trajectories, is ruled out for this proposed RNN, showing that it is stable with respect to changes in the input. This largely on account of the fact that the activation function in B is globally bounded.
18

---

Published as a conference paper at ICLR 2021
E.3 PROOF OF PROPOSITION 3.2
From 6 we readily calculate that,
EE, 0x,
=ly,. -Vn, 0].
(27)
Similarly from 3 we calculate,
A1 Zid5
(A (AR-1)yk-))
At
[A (AR-1)yk-
if (i.j) th entry of W,
1+At m..
mi, mn
A12 Zid (Ax (Ak-1 [k
Zi.
[ k-
)ak
(i. th entry of W,
0 Xx
1+47 mmm
1+227 17.
00
1+At m.d
A12 Zid
(Aa-1)
(
Ar 1+
mdd
if (i.j) th entry of V.
A12 Z
(Ax-1)
At (+AT
Z
i-th entry of b
14At
(28)
where Zm (A+-1) E ma is matrix with all elements are zero except for the (i,j)-th entry
which is o' d' i.e the i-th entry of o' (AA We easily see that |Zi (Ax-1)llo S for all i, J, m, i and all choices of A
Now, using definitions matrix and vector norms and applying 28
10P. together with 27 and
we obtain the following estimate on the norm:
(ly oc + [/.llss +1+ 3(n ")AA"))Ad]]y: bo,
is entry of W,
(
(Ilynll + [[ 3(n k))t") (lAtt))sAr 5 oloo.
if is entry of W,
(29)
00
(lly + [[Tx//==) 3(n k)At") " A uel
if is entry of V, if
(llyxllee + +1151l==) 3(n ())tt")6
is entry of b.
We will estimate the above term, just for the case of 0 is an entry of W, the rest of the terms are very similar estimate.
For simplicity of notation, we let is 1 z k and aim to estimate the term,
< [| + 3(n / - 12 I
3(n A)A
09
I mmVnkAt(11 + 3(n 3(n -0)--)))"11 0A I1 vmw./ + 3(n k) A (by 5
S < + nmv k) + |/ mkv A1 3(n k) At &tt..
(30)
To further analyze the above estimate, recall that nAI In < and consider regimes. Let us start by considering short-term dependencies by letting k = n. i. - k = with constant c (1), independent of n, k. In this case, straightforward application of the above assumptions in the bound (30p yields,
different
de(2)
< mmVnkō + 3m/ +) + | + I
09
< mt,, met S I v
< mo,t + Vm (for At << as a 1/2)
< môAt + Iy
(31)
19

---

Published as a conference paper at ICLR 2021
Next, we consider long-term dependencies by setting &s << n and estimating.
deh!)
I mvnkkss + 3m 3 k) + 18 mm + 33 Vinds
00
S E + 3mt oot + V + 3 mm m
S m 3m5 11 ll má5 + + mm (as <))
s 3m 3 lo Vm (as S 1 and At << 1).
(32)
Thus, in all cases, we have that,
S 35 (m + vmllFnllas) (as r 2 1/2).
(33)
00
Applying the above estimate in 10) allows us bound the gradient by,
iEn
ae
S 3ötn (m + /mllya loo)
(34)
00
00
1Sk5n
Therefore, the gradient of the loss function
6 can be bounded as,
N
00
N
00
n=1
mAt
Am
36
yn 12
N
mAt
m At
38
n
N
N
(35)
3
8(N + A) (m Y/m)
2
8(4N + 4t)(m (
2
3
8(1 Ar) at) + y Vmm
(as =)
2
3
(m 8 8/m),
2
which is the desired estimate 9
E.4 ON THE ASSUMPTION
8
AND TRAINING
Note that all the estimates were based on the fact that we were able to choose a time step At in 3 that enforces the condition (8.. For any fixed weights W, W, we can indeed choose such value of
to satisfy (8). However, we train the RNN to find the weights that minimize the loss function (6).
Can we find a hyperparameter At such that is satisfied at every step of the stochastic gradient descent method for training?
To investigate this issue, we consider simple gradient descent method of the form:
3 Q = 0 - ae
(0x).
(36)
Note that Ç is the constant (non moon-ddapted) learning rate. We assume for simplicity that 6n = 0 (other choices lead to the addition of constant). Then, straightforward estimate on the weight given by.
20

---

Published as a conference paper at ICLR 2021
1et1l
10cl +
00 (0z)
$ 18.1 + 5 (m + Y mvm)
(by (35
(37)
3 1001 + Ec (m Y //m) (c
3 (m + Y F/m).
In order to calculate the minimum number of steps L in the gradient descent method (36) such that the condition (8) is satisfied, we set & = L in (37) and applying it to the condition 8 leads to the straightforward estimate,
2
(38)
[ Y mAI"~58
Note that the parameter å < 1, while in general, the learning rate Ç <<< 1. Thus, as long r S 1, we see that the assumption (8) holds for a large number of steps of the gradient descent method. We remark that the above estimate (38) is large underestimate on L. In the experiments presented in this article, we are able to take a very large number of training steps, while the gradients remain within range (see Fig. 3)
E.5 PROOF PROPOSITION 3
We start with the following decomposition of the recurrent matrices:
Ox;
-M-1 + AtM;-
Ox,-1
A
B 0 0
M,-1 B
M+11=
C,-1
with B, C defined in (2). By the assumption 8 one can readily check that ||M S Al, for
all ls < Sign 1.
We will use an induction argument to show the following representation formula for the product of Jacobians,
Σ At
C
Ox,
0x;
II k<isa
j=k &
+ O(At). (39)
0x
8X,-1
j+1
B
+ Σ
II
C
B
II
C
jmn-2
iin-1
imn-1
We start by the outermost product and calculate,
ax,
0Xn-1 (Mn-1+ ALM (M--2+ AlMn-2)
Ox 0x,22 = -M + A1(M + Mn Mn-1 +AAP))
By direct multiplication, we obtain,
T
At aC 4 C
Mn-1M
B + Cn-IBB-
Cn-1 Cn-2
C,
Bn-
0
At
B Cn-2
Using the definitions in 12 and
we can easily see that
Cn-1BB-2
= -((A).
0
Cn-1
Similarly, it is easy to show that
M, Mn Mn-1 O(A))
21

---

Published as a conference paper at ICLR 2021
Plugging all the above estimates yields,
aXn Ox
At n Cn-) + +0(At))
ax 2xn-2 Bn- Cn
Cn-1
which is exactly the form of the leading term (39
Iterating the above calculations (n k) times and realizing that (n k) nAt At yields the formula (39)
Recall that we have set Wij for some i.j I m in proposition 3.3 Directly calculating
with 27) (28) and the representation formula 39 yields the formula,
ae
00 yn At
SZ m.m (At- y At 8C D mn
+ + O(At")
(40)
with matrix C* defined as,
n-1 k C²
E C J=
and Zini (A Rmxm is matrix with all elements except for the (i, entry which to of (a)-1), i.c. the i-th entry of o' (At-1)..
Note that the formula (40) can be explicitly written as,
8Aioo of +512
8A0 (a oo Σ + O(A(²)
(41)
a0
with y denoting the j-th element of vector and
wayk + Σ W
(42)
f=1
By the assumption (8. can readily see that
|| S1+ At.
Therefore by the fact that o' = sech?, the assumption y = O(vtx) and 422 we obtain,
غ sech 2 ( (VkAt(1+ + At) s da (aa-11 s 1.
(43)
Using 43 in we obtain,
s o (at- : (
(44)
Using the definition of C,, we can expand the product in C' and neglect terms of order O(Ar*), to obtain
b - ( 1)8A7))L
i=j
Summing over j and using the fact that <<< n,, we obtain that
C = (O(n))
(45)
Plugging 45 and 43 into
leads to,
sAt o' (af Σ C = (000211)))
(46)
/=1
Combining 44 and (46) yields the desired estimate (66
22

---

Published as a conference paper at ICLR 2021
Remark. A careful examination of the above proof reveals that the constants hidden in the prefactors of the leading term O (cOAt) of (T6) stem from the formula (46). Here, we have used the
assumption that y' O(vix: Note that this assumption implicitly assumes that the energy bound (5) is equidistributed among all the elements of the vector y& and results in the obfuscation of the constants in the leading term of (16). Given that the energy bound (5) is too coarse to allow for precise upper and lower bounds on each individual element of the hidden state vector yı. we do not any other way of, in general, determining the distribution of energy among individual entries of the hidden state vector. Thus, assuming cquidistribution seems reasonable. On the other hand, in practice, has access to all the terms in formula (46) for each numerical experiment and if one is interested, then one can directly evaluate the precise bound on the leading term of the formula (16).
RIGOROUS ESTIMATES FOR THE RNN 3 WITH ñ = n - 1 AND GENERAL VALUES OF 7
In this section, we will provide rigorous estimates, similar to that of propositions 3.1 in (3) leading to,
E.1 and 3.2 for
the version of coRNN (3 that results by setting
Yn Yn-1 + Atzn,
(47)
n 2n-1 Ato (Wyn-1 Wa. Vu, b) Atyyn-1 Ateza-1-
Note that
(47 can be equivalently written as,
yn -n Alzn,
(48)
Zn = (1 <Lt) 2n-1 Ato (Wyn-1 WZ.-1 Vu, + +b) - -ttyn-1-
We will also consider the case of non-unit values of the control parameters 7 and below.
Bounds on Hidden states. We start the following bound on the hidden states of 47).
Proposition F11 Let the damping parameter e > following condition,
2 and the time step At in the RNN sctisfy the
2c
At
(49)
7 + 2
Let yn; 2n, be the hidden states of the RNN following (energy) bounds:
47 for
S S N, then the hidden states satisfy the
1 YnYn -2, Zn
mta 7
(50)
We set A An-1 Wyn-1 + Wz Vu,- + and as in the proof of proposition 3.1 we multiply
(y , zz, to (47) and use elementary identities and rearrange terms to obtain,
y Yn
2
Nn-1
2T 1hn-
(yn yn -n-1) (yn yn-1)
27
27
2
(zn Zn (a-1))(zn 2nn11)
27
At 2,2 (o(Ann))))
cAt
eAt 2n (2n - 2 2n-1)-
7 2, Un
We use rescaled version of the well-known Cauchy's inequality
ca² 2 b²
ab Σ 2
2c
for constant to be determined, rewrite the above identity as,
vTTY Yn
T
1--1n
2 1n-1 + (yn. yn (yn yn-1)
2y
24
CAt
1
At
(Zn 2n-1) (zn, T
Zn-2) +
o(A (An-1)) To (An-1)
207
27
24
AI
Ce Al.
cAl
T
2, Znn
24
2-
7
23

---

Published as a conference paper at ICLR 2021
Using the first equation in
(47
the above inequality reduces to,
Yn yIyy Y
2n
18n-1
Pn-1
27
24
cAt
1
At
(Zn 2n-1) (z, kn-1) + o(nn-)) (Ann1) 24
207
27
At
At
CE
eAt
T
+
3n Zn.
2
2x
2
7
As long as,
c (2 - c)c - 1)
At S min
(51)
c'
?
we can easily check that,
Yn yn
J. Zn
Yn 1n-1
Z 2n-1Zn-1 -n-1
At
o ((An o((n-1)) (Ao--))T((Annii))
2
2%
2
2y
2y
Yn Yn-1In-
1- 2 n-
An-14n-1
1n-1
mAt
(o 1).
2
2
2
Iterating the above bound till n = and using the zero initial data yields the desired (50) as long as we find such that the condition (5T) is satisfied. To do we equalize the two terms on the right hand side of 51 to obtain,
<(2c 1)
? 7 E2
From the assumption 49) and the fact that 2. value of >> and (51 is satisfied, which completes the proof.
that such > 0 always exists for any
We remark that the same bound on the hidden states is obtained for both versions of coRNN, i.e. 3 with î = n and (47). However, the difference lies in the constraint on the time step A.. In contrast to (49). a careful examination of the proof of proposition 3.1 reveals that the condition on the time step for the stability of (3) with π = n is given by.
At 2E - 1
(52)
and is clearly less stringent than the condition 5T) for the stability of (47). For instance, in the prototypical case y = the other hand, the stability of (47) is ensured as long as At < 2. However, it is essential to recall that these conditions are only sufficient to ensure stability and are by no means necessary. Thus in practice, the coRNN version 447 is found to be stable in the same range of time steps as the version 33 with n.
the stability of (3) with
ensured for any At
On
On the exploding and vanishing gradient problems for CoRNN 47 Next, we have the following upper bound on the hidden state gradients for the version (47 of CORNN,
Proposition F.2 Let yn Zn be the hidden states generated by the RNN (47). We assume that the damping parameter e > 2 satisfies,
I and the time step At can be chosen such that in addition to (51) it also
max (At(y + [[]])s) = ČAF
11sl,
(53)
and with the constant Ò independent of the other parameters of the RNN 47 Then the gradient of
the loss function & (6) with respect to any parameter O E + is bounded as,
36
3(Ō) (m + Y Vm)
<
(54)
ae
27
with the constant Č. defined in 53 and Y max data
|| be a bound on the underlying training N
24

---

Published as a conference paper at ICLR 2021
The proof of this proposition is completely analogous to the proof of proposition 3.2 and omit the details here.
Note that the bound (54) enforces that hidden state gradients cannot explode for version (47 of CoRNN. similar statement for the vanishing gradient problem inferred from the proposition below.
Proposition F.3 Let yn be the hidden states generated by the RNN 47
Under the assumption that for all 1 iis m and 53). the gradient for long-term dependencies satisfies,
Y = (V tm 7
( A= +0
è 5(1 5)att +O(AP), ē= sech? (VRAt(1+ + At)
08 = 0
)
k k<<n n.
(55)
The proof is repetition of the steps of the proof of proposition 3.3 with suitable modifications for the structure of the RNN and non-unit E, yy and we omit the tedious calculations here. Note that (55) rules the vanishing gradient problem for the CoRNN version 47).
25


Reviews:
[{"review_id": "t6jo73VdY9Z", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Firstly, this paper conducts the rigorous analysis of the coRNN via the formula deduction to verify the bound. Then the coRNN is proved to mitigate the exploding and vanishing gradient problem and this is also validated in a series of experiments. Also, the performance of the coRNN is comparable or better compared to state-of-the-art models. This paper provides a new idea to address the exploding and vanishing gradient problem, which hinders the development of deeper neural networks tremendously. In my opinion, this coRNN model is meaningful for practical application, especially for the extension of more complicated neural networks. \n\nBesides, for the biomedical signals with high temporal resolution (e.g., electroencephalogram, electromyogram), the coRNN model can be a good alternative in future work. Furthermore, the efficiency of the proposed model also be proved by the mathematic formulation and experiments ranging from pure synthetic tasks designed to learn long-term dependencies to more realistic tasks rigorously. Considering the whole structure of this paper, I argue that the clarity is clear and logical. Different from the recently published literature, this paper has explicit use of networks of oscillators with the underlying biological motivation, so this paper expresses the originality in some extent. To sum up, the quality of this paper is suitable for the publication in ICLR2021. \n\nThere are two main pros in this paper: \n1.\tThe theoretical verification is clear and rigorous, readers can easily catch good understanding of the bounds this paper proves following the formula deduction. Specifically, this paper demonstrates how to avoid the exploding and vanishing gradient problem for the RNN in theory.  \n2.\tThe experiments are quite abundant, experimental results show that the coRNN can not only avoid the exploding and vanishing gradient problem, but also achieve better performance with fewer parameters compared to recent studies.  \n\nBut some cons should also be noticed. Firstly, the illustration of proposed coRNN should be presented in the paper, which is more comprehensible. Secondly, the related work part, when mentioning the similar works, it will be better to describe the main differences and correction with this paper more specifically. Lastly, in the part of discussion, the practical significance of proposed coRNN should be emphasized with more words.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "ngsOAO-EcZM", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes a novel RNN architecture (CorNN) to tackle the infamous problem of vanishing and exploding gradients in RNNs. The novel CorNN architecture is based on time-discretized forced coupled damped nonlinear oscillators. For the gradient norm of CorNN analytical lower and upper bounds are calculated implying that CorNN avoids vanishing and exploding gradients.\nThis is accompanied by numerical results (including code) that demonstrate the improved trainability of CorNN in permuted sequential MNIST, and adding task and noise-padded CIFAR10 compared to some other RNN architectures (GRU, LSTM, antisymmetricRNN, IMDB sentiment analysis, and a human activity recognition task.\n\nIn summary, the paper proposes a useful and mathematically transparent way of tackling the challenge of training RNN on tasks with long-time dependencies. \n\nWeak points of the paper:\n* In the mathematical analysis the paper claims to \"rigorously prove precise bounds on [...] gradients, enabling the solution of the exploding and vanishing gradient problem\". If I am not misunderstanding, the mathematical bounds are actually only shown for the initial gradient norm (plus some number of steps afterward in section C4 of the Supplementary Material).\n\n* A more in-depth comparison of numerical gradient norms and their respective upper/lower analytical bounds would be desirable, as the lower bound on the gradient in (16) is only given in as O(∆t^(3/2)) without any prefactors. \n\n* While the numerical results demonstrate improved trainability and performance on a number of tasks, the underlying reasons remain mostly unclear. Supplementary Material B gives some heuristics on how a superposition of forced coupled damped nonlinear oscillator can generate complex output, but it doesn't explain the avoidance of exploding and vanishing gradients and the superiority to other solutions (e.g. gated units like LSTM or GRU).\n\n* Missing: What are the limitations of CorNN, when would you expect them to fail?\n\n* The link to biological networks is not yet very convincing. While there exist without doubt many oscillations on different scales in the brain, it is not clear how the insights gained here could be applied to the brain.\n\nSome smaller comments:\n* How could expressivity be quantified?\n* figure 1: Plotting MSE with a logarithmic or semilogarithmic axis would help to distinguish small errors.\n* figure 3: It would help to also plot this for other tasks.\n* figure 3 lines for other tasks would be helpful", "extra_scores": {}, "extra_text": {}}, {"review_id": "it96G_j8Fwv", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a new continuous-time formulation for modeling recurrent units. The particular form of the recurrent unit is motivated by a system of coupled oscillators. These systems are well studied and widely used in the physical, engineering and biological sciences. Establishing this connection has the potential to motivate interesting future works. The performance of the proposed recurrent unit is state of the art.\n\n\nReasons for my score: Overall, I vote for marginally above acceptance threshold. I like very much the proposed approach for modeling recurrent units. Further, the presented results are intriguing, and the paper is well written. However, I have some concerns (see below). I am happy to increase my score if the authors can address my concerns in the rebuttal period. \n\n\nPros: \n-------\n+ Second-order systems of ODEs seem to be a promising approach for modeling recurrent units, and this approach has not received much attention for this task before. Indeed, this paper impressively demonstrates that a unit motivated by a system of coupled oscillators is able to achieve state of the art performance on a range of benchmark tasks.\n\n+ The analysis shows that the particular form of the proposed continuous-time unit mitigates the vanishing and exploding gradients problem by design, which is very appealing. The analysis is mathematically sound.\n\n+ Code is provided!\n\n\nCons: \n-------\n- In Eq. (3) the authors advocate the IMEX scheme to obtain a discretization of (2). I was very curious to see how the authors implement this scheme in practice, however, the provided implementation revealed that the authors use an explicit scheme in practice. Please, comment why you chose the IMEX scheme here. Does the analysis also hold if you use an explicit discretization in (3), and if, why do you mask the fact that you are using an explicit scheme in practice. I feel, it would be relevant to discuss how you train the unit in practice.   \n\n- Section 3 is no pleasure to read. It is not clear to me what the value of the sketch of the proofs are, since the proofs are pretty standard. Instead, the space could be better used for an extended qualitative discussion of the analysis and the nice properties of the proposed recurrent unit. For instance, you can extend the discussion around proposition 3.1 and provide some context on why you want to rule out chaotic behavior (this might not be obvious for everyone); further it would be nice to see a better discussion on the effect of dampening and forcing on the performance of the recurrent unit. Also, I would like to suggest to move parts of Appendix B into the main text, since this discussion actually helps to build some intuition for the proposed unit.\n\n- The extremely good performance on the sMNIST task is slightly surprising, since my intuition would not suggest that the particular form of the unit has the ability to substantially improve the expressivity as compared to some other recently proposed units. I used the provided code to evaluate the coRNN (N=256 and 128) on the sMNIST task and the highest accuracy that I was able to obtain (out of 8 runs on 4 different GPUs) was 99.2% on the test set. These results are still very good, but they do not match the reported results. (Note, that the code is printing out the accuracy for a smaller validation set which indicates a higher accuracy than is actually obtained on the test set.) This said, I would like to ask the authors to double check the experiments on sMNIST. (Also, I assume that the model can be trained in less time if the learning rate is decayed much earlier, e.g., around epoch 30 and a second time around epoch 60.)\n\n- It is not clear to me how sensitive the RNN is to the particular choices of \\gamma and \\epsilon. It would be good to provide some form of ablation study that studies how the performance varies for different values of \\gamma (and \\epsilon) while keeping all other tuning parameters fixed (I assume that you have all these results handy since you have performed an extensive hyperparamter search). This would help to gain some better intuition for how difficult it is to tune the proposed unit. In other words, I would like to see how sharp the performance drop is if you perturb the tuning parameters slightly (i.e., plot the test accuracy as a function of \\gamma and \\epsilon).\n\n\nMinor comments: \n-------\n* Given additional space, it would be nice to see an extended related work section.\n\n* It would be nice to so results for a language modeling task.\n\n* In Table 2, the citation for the Fast RNN is incorrect.", "extra_scores": {}, "extra_text": {}}, {"review_id": "NbD6KYwtnPY", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper introduces a novel recurrent neural network architecture which approximately preserves the norm of the gradient irrespective of the number of unroll steps. This complements a rapidly growing line of research that aims to better understand dynamical properties of RNNs and their gradients, thus potentially enabling training of models that capture long term dependencies while avoiding exploding and vanishing gradients.\n\nThe submission has the following strengths to it:\n\n1. It offers a clear and succinct proof of the gradient stability as well as the stability of the forward dynamics.\n1. It provides convincing experiments on relevant datasets, all while showing competitive results.\n1. It's exceedingly well written, and readily understandable even without prior knowledge.\n\n##### I firmly believe that based on those, it should be accepted. It has all the hallmarks of a good, paper with potentially wide application.\n\nThat being said, I do have some minor objections:\n\n1. When presenting Eq. (1), the 'intuitive' interpretation of $\\gamma, \\epsilon$ should given right away, rather than deferred to later sections of the paper, especially the appendix.\n1. The motivation for the use of non-linear oscillators is well-written but perhaps should be de-emphasized. I would like to put-forward the following argument for it. It appears that the choice of the dynamics only constitutes half of the 'puzzle'. \nThe choice of IMEX is mentioned *en passant*, but seems to be rather crucial to obtaining the theoretical guarantees viz. Proposition 3.2 and 3.3. I did not have time to check the derivation against other schemes, but I presume that the choice of IMEX was highly non-trivial in designing the new architecture. If that is indeed the case, then the role of the solver ought to be emphasized.\n1. I appreciate the authors comments regarding the expressivity of the proposed architecture, as well as the demonstrations in Appendix B.\nHowever, I would also appreciate a simple example of the kind of dynamics where the coRNN cell ought to break down -- given the corollary of Proposition 3.1, it would be interesting to show the potential break-down of the network on task that involves approximating a chaotic dynamical systems. In particular, it would be very interesting how coRNN fares against similarly sized gated RNN (LSTM or GRU).\n1. For my own understanding, is the exponent $r$ in Eq. 8 there only to conveniently related $\\eta$ to $\\Delta t$? If not, does it admit some more intuitive interpretation?\n1. Since, the proposed architecture requires a sufficiently small step-size, are the resultant equivalent (in some suitable sense be it topologically, having the same invariant set) to the continuous time dynamics?\n1. Lastly, the authors report the hidden unit dimensionality, but from the main text it is entire unclear whether that's the dimensionality of $y$ or $ \\left[y^\\top, z^\\top \\right]^\\top$. Having looked at the code, it appears to be the former.\n\n#### Edit:\nOut of curiosity I ran the submitted code for the permuted sequential MNIST task, and noticed that the following:\n1. The numbers that the authors report in the paper seem to be result of a single network realization. While granted, this is somewhat consistent with practices common in the community, it makes one question how representative they are of different initial seeds. In the current setting it's hard to disambiguate whether the random seed was chosen coincidentally or rather specifically because of the purported state-of-the-art outcome.\n1. For this reason I suggest the authors compute additional iterates of the model and report some distributional information about the best loss/accuracy for all the tasks covered in the submission.\n1. Attaining \"state-of-the-art\" results is notable, but is by no means pre-requisite for this to be considered a good submission.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: A novel second order nonlinear oscillator RNN architecture is proposed, analyzed, and evaluated in this paper. The results are solid and impactful. Authors and expert reviewers showed exemplary interactions with each other, improving the manuscript in significant ways. All four reviewers overwhelmingly recommended accept. I recommend that this paper be selected as an oral presentation.
Decision: Accept (Oral)