# Recurrent Neural Networks for Dynamical Systems: Applications to Ordinary
Differential Equations, Collective Motion, and Hydrological Modeling

Y. Park  K. Gajamannage<EMAIL_ADDRESS>Department of
Mathematics & Statistics, Texas A&M University – Corpus Christi, Corpus
Christi, TX 78412, USA  D.I. Jayathilake  Department of Physical &
Environmental Sciences, Texas A&M University – Corpus Christi, Corpus Christi,
TX 78412, USA  E. M. Bollt  Department of Electrical and Computer Engineering,
Clarkson University, Potsdam, NY 13699, USA

(August 27, 2024)

###### Abstract

Classical methods of solving spatiotemporal dynamical systems include
statistical approaches such as autoregressive integrated moving average, which
assume linear and stationary relationships between systems’ previous outputs.
Development and implementation of linear methods are relatively simple, but
they often do not capture non-linear relationships in the data. Thus,
artificial neural networks (ANNs) are receiving attention from researchers in
analyzing and forecasting dynamical systems. Recurrent neural networks (RNN),
derived from feed-forward ANNs, use internal memory to process variable-length
sequences of inputs. This allows RNNs to applicable for finding solutions for
a vast variety of problems in spatiotemporal dynamical systems. Thus, in this
paper, we utilize RNNs to treat some specific issues associated with dynamical
systems. Specifically, we analyze the performance of RNNs applied to three
tasks: reconstruction of correct Lorenz solutions for a system with a
formulation error, reconstruction of corrupted collective motion trajectories,
and forecasting of streamflow time series possessing spikes, representing
three fields, namely, ordinary differential equations, collective motion, and
hydrological modeling, respectively. We train and test RNNs uniquely in each
task to demonstrate the broad applicability of RNNs in reconstruction and
forecasting the dynamics of dynamical systems.

††preprint: AIP/123-QED

> Learning non-linear systems in nature requires researchers to solve
> complicated mathematical equations which often involve both high
> computational complexity and less computational precision [1]. Artificial
> neural networks (ANNs) are a class of algorithms that model underlying
> relationships of a set of data through a process that mimics the way the
> human brain operates. Data-driven approaches such as ANNs are capable of
> affording such nonlinear systems with both high efficiency and good
> accuracy. ANN models are increasingly popular for two reasons. First, they
> guarantee that arbitrary continuous dynamical systems can be approximated by
> ANNs with a sufficient amount of sub-models which are called hidden layers
> [2]. Second, ANNs themselves can be viewed as discretizations of continuous
> dynamical systems which make them suitable for studying dynamics [3]. A
> recurrent neural network (RNN) is a class of ANNs where the connections
> between nodes form a directed graph along a temporal sequence where the
> output of a hidden node at the current time-step is sent to the
> corresponding hidden node for the next time-step. Among all the ANN
> architectures, RNNs are considered to be the most faithful candidates for
> modeling and forecasting temporal data sequences such as time series data.
> Thus, here we implement RNN to reconstruct and forecast data generated by
> several nonlinear dynamical systems.

##  I Introduction

Dynamical systems may exhibit dynamics that are highly sensitive to initial
conditions, which is referred to as the butterfly effect. As a result of this
sensitivity, which often manifests itself as an exponential growth of
perturbations in the initial conditions, the behavior of dynamical systems is
sophisticated [4]. From a historical point of view, reconstruction or forecast
of dynamical systems, especially those with spatiotemporal solutions, has been
based on the techniques like linear parametric autoregressive models, moving-
average models, autoregressive moving-average models, etc. All these models
are linear and some of them possess the assumption that the data is being
stationary, thus, these models are not able to cope with complex relationships
in the data [5].

We utilize RNNs with unique ways of training and testing routings for modeling
a general class of dynamical systems since they are considered to be the
state-of-the-art machine learning methods for sequential data. Applications of
RNNs include image classification, object recognition, object detection,
speech recognition, language translation, voice synthesis, etc [6]. An ANN
usually involves a large number of parallelly operating processing units that
are arranged in layers. The first layer receives the raw input information
analogous to optic nerves in human visual processing. Each successive layer
receives the output from both the layer preceding it and the same layer,
rather than from the raw input, as analogous to the way that the neurons
further from the optic nerve receive signals from those closer to it [6]. The
last layer produces the output of the system. RNN possesses an internal state
or short-term memory due to the recurrent feedback connections, as seen by
directed loops in Fig. 1. This makes RNN is suitable for modeling sequential
data, especially time series. Moreover, these connections allow well-trained
RNNs to regenerate dynamics of any temporal nonlinear dynamical system up to a
satisfactory accuracy. Thus, RNN models are widely used to analyze diverse
time series problems including dynamical systems [7]. Most of the complex RNN
architectures, such as Long Short-Term Memory (LSTM) [8] and Gated Recurrent
Unit [9], can be interpreted as a variation or as an extension of the basic
RNN scheme[10].

Figure 1: Standard single-layer recurrent neural network (RNN), i.e., one-
stacked RNN, architecture. The right-hand side schematic is the unrolled
version, with respect to time, of the left-hand side RNN architecture where
the directed loop represents recurrent feedback between different states of
the same layer. Here, $X$, $\boldsymbol{w}_{hh}$, and $\hat{Y}$ represent
input, the weights of the recurrent connections, and the output. Moreover, for
the $t$-th state, while $\boldsymbol{h}^{(t)}$ represents the hidden state,
$w^{(t)}_{ih}$, $w^{(t)}_{hh}$, and $w^{(t)}_{ho}$ represent weights between
input and hidden layer, states of the hidden layer, and hidden layer and
output, respectively. Note that, each state itself is a conventional
artificial neural network; however, the RNN is formed by linking consecutive
states within each layer with a directed edge. This recurrent feedback accords
memory for an RNN through its states (i.e., time-steps).

The objective of this study is to demonstrate how well RNNs learn
spatiotemporal dynamics of multi-scale dynamical systems. Our focus is either
reconstruct or forecast short-term or long-term trajectories of a dynamical
system of interest with given initial conditions. We use three exemplary
datasets that are generated from three dynamical systems, namely, the Lorenz
attractor [11], a generalized Vicsek model [12], and a streamflow model, where
those represent three diverse fields, ordinary differential equations,
collective motion, and hydrological modeling, respectively. For given orbits
that are generated from a Lorenz system with a formulation error, we train an
RNN to eliminate the formulation error and reconstruct the correct system’s
responses. Then, we generate a collective motion dataset using the generalized
version in [12], of the classic Vicsek model in [13], and create two copies of
that, one by imposing some additive Gaussian noise and the other by deleting
some segments of the trajectories. Here, our goals in these experiments are to
eliminate the additive noise from the noisy trajectories. Finally, we forecast
streamflow data that are naturally highly volatile and consist of frequent
spiky fluctuations.

This paper is structured as follows: In Sec. II, we present the architecture
and the optimization routine of RNNs for analyzing and forecasting data
generated from spatiotemporal dynamical systems. In Sec. III, we validate the
performance of RNNs applied to three dynamical system datasets representing
three fields, ordinary differential equations, collective motion, and
hydrological modeling. The analysis conducted using RNNs addresses three
unique issues and solutions for them. Sec. IV states the discussion including
key findings, limitations, future work, and conclusions. Table 1 represents
notations used in this paper along with their descriptions.

Table 1: Notations used in this paper and their descriptions Notation |  Description
---|---
$N_{i}$ |  Number of nodes in the input layer of an RNN.
$N_{h}$ |  Number of nodes in the hidden layer of an RNN.
$N_{o}$ |  Number of nodes in the output layer of an RNN.
$N$ |  Total number of orbits, trajectories, and agents.
$K$ |  Total number of layers.
$t$ |  Index for time-steps where $1\leq t\leq T.$
$k$ |  Index for the RNN’s layers where $1\leq k\leq K.$
$n$ |  Index for orbits, trajectories, or agents where $1\leq n\leq N.$
$\sigma,\rho,\beta$ |  Parameter in Lorenz.
$T$ |  Number of total time-steps
$\delta$ |  Step size of the generalized vicsek model
$\eta$ |  Corruption level of orbits.
$\sigma$ |  Standard deviation of Gaussian noise distribution.
$t$ |  Index for time-steps where $1\leq t\leq T.$
$\epsilon^{(t)}_{i}$ |  Noise parameter imposed to the $i$-th agent at the $t$-th time-step.
$\theta^{(t)}_{i}$ |  Orientation of the $i$-th agent at the $t$-th time-step.
$r_{d}$ |  Radius of interaction
$r^{(t)}$ |  Radius of spiral example.
$k^{(t)}$ |  Angle of rotation of the a spiral at time $t$ from the initial position.
$N^{(t)}_{i}$ |  $r$-radius of neighborhood of the $i$-th agent at the $t$-th time-step.
$\boldsymbol{x}^{(t)}$ |  Input vector at time-step t.
$\boldsymbol{h}_{k}^{(t)}$ |  Hidden state of the $k$-th layer for the $t$-th time-step.
$\boldsymbol{y}^{(t)}$ |  Label vector at time-step t.
$\boldsymbol{b}_{i}$ |  Bias vector for the input layer.
$\boldsymbol{b}_{h}$ |  Bias vector for the hidden layer.
$\boldsymbol{x}^{(t)}_{n}$ |  Input vector of $n$-th orbit at time-step $t$ ($t$-th point of $n$-th orbit).
$\boldsymbol{y}^{(0)}_{n}$ |  Initial condition of $n$-th orbit.
$\boldsymbol{y}^{(t)}_{n}$ |  Label vector of the $n$-th orbit at time-step t.
$\hat{\boldsymbol{y}}^{(t)}_{n}$ |  Output vector of the $n$-th orbit at time $t$ corresponding to $\boldsymbol{x}^{(t)}_{n}$.
$\boldsymbol{u}^{(t)}_{i}$ |  Average direction of motion of the agents in the neighborhood of $N^{(t)}_{i}$.
$\boldsymbol{c}^{(t)}$ |  The $t$-th point on the centroid of the spiral.
$\boldsymbol{v}^{(t)}_{i}$ |  Velocity of the $i$-th agent at the $t$-th time-step.
$X$ |  Set of alll input data.
$\hat{Y}$ |  Set of all output.
$Y$ |  Set of all label.
$f$ |  Linear transformation in RNN.
$g$ |  Non-linear transformation in RNN.
$L(\cdot)$ |  Loss function.
$W_{ih}$ |  Weight matrix between the input layer $i$ and the hidden layer $h$.
$W_{h_{k}h_{k^{\prime}}}$ |  Weight matrix between the hidden layer $h_{k}$ and the hidden layer $h_{k^{\prime}}$.
$W_{ho}$ |  Weight matrix between the hidden layer $h$ and the output layer $o$.
$R^{(t)}_{j}$ |  Rotation matrix of the $j$-th agent at the $t$-th time-step.

##  II Recurrent neural networks

RNNs are a major class of machine learning tools that are particularly well
suited to work on sequential data. Connections between nodes form a directed
graph along a temporal sequence, that can be used to learn nonlinear
dependencies of a system over time [14]. In its most general form, an RNN can
be seen as a weighted directed graph that contains three different kinds of
nodes, namely, the input nodes, hidden nodes, and output nodes [15]. As seen
in Fig. 1, where we represent the structure of a simple RNN with one layer,
input nodes do not have incoming connections and output nodes do not have
outgoing connections, but hidden nodes have both. By RNN’s design, two
different nodes which are either at the same time level or at different time
levels can be connected by an edge, as seen in Fig. 2.

Figure 2: Standard stacked recurrent neural network (RNN) architecture of
having three layers, i.e., three-stacked RNN. The right-hand side schematic is
the unrolled version, with respect to time, of the left-hand side RNN
architecture where directed loops represent recurrent feedbacks between
different states of the same layer. Here, $X$, $\boldsymbol{w}_{h_{k}h_{k}}$,
and $Y$ represent input, the weights of the recurrent connections of the
$k$-th layer, and the output. For the $t$-th state, where $1\leq t\leq T$ for
some $T$, $\boldsymbol{x}^{(t)}$, $\boldsymbol{h}_{k}^{(t)}$, and
$\hat{\boldsymbol{y}}^{(t)}$ represent the input vector, the hidden state of
the $k$-th layer where $1\leq k\leq K$ for some $K$, and the output vector,
respectively. Moreover, $w^{(t)}_{ih_{1}}$, $w^{(t)}_{h_{k}h_{k^{\prime}}}$,
and $w^{(t)}_{h_{3}o}$ are the weights between input and hidden layer 1,
hidden layer $h_{k}$ and hidden layer $h_{k^{\prime}}$, and hidden layer 3 and
output, respectively. For the $k$-th layer, $w^{(t)}_{h_{k}h_{k}}$ represents
the weight between two consecutive states.

An RNN has a memory that stores _knowledge_ about the data sequence that it
has already seen, however, its memory of internal states through iterations is
short-term [16]. By applying activation functions (also called, transfer
functions) to previous states and inputs, RNNs compute new states. Typical
activation functions in RNNs include sigmoid, hyperbolic tangent, and ReLU;
however, sigmoid and hyperbolic tangent encounter vanishing of the gradient in
the optimization step [17]. Properly customized RNNs with the best activation
function and the best number of hidden nodes can capture the dynamic of any
nonlinear system up to a high precision [18]. The model capacity and
flexibility for learning non-linear systems can be increased by a multi-layer
recurrent structure [19]. By this multi-layer structure, the hidden output of
one recurrent layer can be propagated through time and is used as the input
data to the next recurrent layer. In our study, we used three-stacked RNNs (or
recurrent layers) as seen in Fig. 2.

###  II.1 Architecture of RNNs

For simplicity, we formulate an RNN with only one hidden layer, i.e., one-
stacked RNN, as shown in Fig. 1. We denote the number of nodes in the input
layer, hidden layer, and output layer as $N_{i}$, $N_{h}$, and $N_{o}$. For
the $t$-th state, the weight matrix from input to hidden layer is denoted as
$w^{(t)}_{ih}\in\mathbb{R}^{N_{i}\times N_{h}}$, the weight matrix between
hidden layer as $w^{(t)}_{hh}\in\mathbb{R}^{N_{h}\times N_{h}}$, and the
weight matrix from the hidden layer to the output layer as
$w^{(t)}_{ho}\in\mathbb{R}^{N_{h}\times N_{o}}$. Moreover, for the $t$-th
state, we denote the bias vector for the input layer as
$\boldsymbol{b}^{(t)}_{i}\in\mathbb{R}^{N_{i}}$ and the bias vector for the
hidden layer as $\boldsymbol{b}^{(t)}_{h}\in\mathbb{R}^{N_{h}}$. Here, we
develop the formulation for the $t$-th state, where $1\leq t\leq T$ for some
$T$, i.e. for the $t$-th ANN, in the RNN. The input $\boldsymbol{x}^{(t)}$ of
the $t$-th state is multiplied with $w^{(t)}_{ih}$, and summed that with
$\boldsymbol{b}^{(t)}_{i}$, i.e.
$w^{(t)}_{ih}\boldsymbol{x}^{(t)}+\boldsymbol{b}^{(t)}_{i}$, to get the flow
of information from the input to the hidden layer. Similarly, the information
flow by the recurrent feedback from the hidden node of the ($t-1$)-th state to
that of the $t$-th state is
$w^{(t-1)}_{hh}\boldsymbol{h}^{(t-1)}+\boldsymbol{b}^{(t-1)}_{h}$. The hidden
state $\boldsymbol{h}^{(t)}$ of the $t$-th state is the sum of the above two
information flows passing through a nonlinear transformation, say $f$, such
that,

$\boldsymbol{h}^{(t)}=f\left(w^{(t)}_{ih}\boldsymbol{x}^{(t)}+\boldsymbol{b}^{(t)}_{i}+w^{(t-1)}_{hh}\boldsymbol{h}^{(t-1)}+\boldsymbol{b}^{(t-1)}_{h}\right),$
(1)

where $f$ is the activation function of the neurons which is usually sigmoid
or hyperbolic tangent. The information flow from the hidden layer of the
$t$-th state to its output, i.e.
$w^{(t)}_{ho}\boldsymbol{h}^{(t)}+\boldsymbol{b}_{o}$, is composed with
another transformation, say $g$, to get the output, say
$\hat{\boldsymbol{y}}^{(t)}$, of the RNN such that,

$\hat{\boldsymbol{y}}^{(t)}=g\left(w^{(t)}_{ho}\boldsymbol{h}^{(t)}+\boldsymbol{b}^{(t)}_{o}\right),$
(2)

where $g$ is linear in general. For all $t$’s, $\boldsymbol{h}^{(t)}$ is
usually initialized with a vector of zeros. The explicit correspondence of
Eqns. (1) and (2) with the RNN’s architecture, helps us adopt spatiotemporal
notion of RNNs using so-called _finite unfolding in time_ with shared weight
matrices [20].

By learning the formulation pattern in Eqns. (1) and (2), we derive the
formulation for a stacked RNN with $K$ hidden layers. Note that, Fig. 2
provides the architecture for such an RNN with $K=3$, i.e., three-stacked RNN.
For the $t$-th time-step, let $\boldsymbol{x}^{(t)}$,
$\boldsymbol{h}_{k}^{(t)}$, and $\boldsymbol{y}^{(t)}$ represent the input
vector, the hidden state of the $k$-th layer, and the output vector,
respectively. We denote the number of nodes in the input layer, $k$-th hidden
layer, and output layer as $N_{i}$, $N_{h_{k}}$, and $N_{o}$. Moreover,
$w^{(t)}_{ih_{1}}\in\mathbb{R}^{N_{i}\times N_{h_{1}}}$,
$w^{(t)}_{h_{k}h_{k^{\prime}}}\in\mathbb{R}^{N_{h_{k}}\times
N_{h_{k}^{\prime}}}$, and $w^{(t)}_{h_{3}o}\in\mathbb{R}^{N_{h_{3}}\times
N_{o}}$ are the weights between the input and hidden layer 1, hidden layer
$h_{k}$ and hidden layer $h_{k^{\prime}}$, and hidden layer 3 and the output,
respectively. Similar to Eqn. (1), the $t$-th state of the first hidden layer
of the RNN is sent through a nonlinear activation function $f$ along with both
the input $\boldsymbol{x}^{(t)}$ and the $(t-1)$-th state of the first hidden
layer, $\boldsymbol{h}^{(t-1)}_{1}$, as

$\boldsymbol{h}^{(t)}_{1}=f\left(w^{(t)}_{ih_{1}}\boldsymbol{x}^{(t)}+\boldsymbol{b}^{(t)}_{i}+w^{(t-1)}_{h_{1}h_{1}}\boldsymbol{h}_{1}^{(t-1)}+\boldsymbol{b}^{(t-1)}_{h_{1}}\right),$
(3)

where $\boldsymbol{b}^{(t)}_{i}$ is the bias vector for the input at the
$t$-th state and $\boldsymbol{b}^{(t-1)}_{h_{1}}$ is the bias vector for the
first hidden layer at the $t$-th state. The hidden state,
$\boldsymbol{h}^{(t)}_{k}$ for $k=2,3,\dots,K$ of the $t$-th state is based on
both its $t$-th state at the previous hidden layer, i.e.,
$\boldsymbol{h}^{(t)}_{k-1}$, and the $(t-1)$-th state of the same hidden
layer, i.e., $\boldsymbol{h}^{(t-1)}_{k}$, as

$\boldsymbol{h}^{(t)}_{k}=f\left(w^{(t)}_{h_{k-1}h_{k}}\boldsymbol{h}^{(t)}_{k-1}+\boldsymbol{b}^{(t)}_{h_{k-1}h_{k}}+w^{(t-1)}_{h_{k}h_{k}}\boldsymbol{h}^{(t-1)}_{k}+\boldsymbol{b}^{(t-1)}_{h_{k}h_{k}}\right),$
(4)

where $\boldsymbol{b}^{(t)}_{h_{k-1}h_{k}}$ is the bias vector for the
$(k-1)$-th hidden later at the $t$-th state and
$\boldsymbol{b}^{(t-1)}_{h_{k}h_{k}}$ is the bias vector for the $k$-th hidden
layer at the $(t-1)$-th state. The output of the $t$-th state, say
$\hat{\boldsymbol{y}}^{(t)}$, is only based on its $K$-th hidden layer such
that,

$\hat{\boldsymbol{y}}^{(t)}=g\left(w^{(t)}_{h_{K}o}\boldsymbol{h}^{(t)}_{K}+\boldsymbol{b}^{(t)}_{h_{K}o}\right),$
(5)

where $\boldsymbol{b}^{(t)}_{h_{K}o}$ is the bias vector for the $K$-th hidden
layer at the $t$-th state. The stacked RNN’s initial hidden states
$\boldsymbol{h}^{(0)}_{k},\dots,\boldsymbol{h}^{(T)}_{k}$ should be
initialized for each layer $k$. All the weight matrices and bias vectors are
optimized using the gradient descent method according to a straightforward
backpropagation through time (BTT) procedure [21].

###  II.2 Optimization of RNNs

Training an RNN is mathematically implemented as a minimization of a relevant
reconstruction error, widely called as the loss, function with respect to
weights and bias vectors of Eqns. (3), (4), and (5). This optimization is
carried out in four steps: first, forward propagation of input data through
the neural network to get the output; second, calculate the loss between
forecasted output and the expected output; third, calculate the derivatives of
the loss function with respect to the ANN’s weights and bias vectors using
BTT; and fourth, adjusting the weights and bias vectors by gradient descent
method [22].

BTT unrolls backward all the dependencies of the output on the weights of the
system [23], which is represented from left side to right side in Figs. 1 and
2. We train the RNN by one data instance, say instant $n$, at a time; thus, we
formulate the loss function of RNN’s as a normalization over all the data
instances, say $N$. We denote the $n$-th input data instance as
$\boldsymbol{x}_{n}=[\boldsymbol{x}^{(1)}_{n},\dots,\boldsymbol{x}^{(T)}_{n}]\in
X$ where $X$ represents the set of all the input data, and the corresponding
output as
$\hat{\boldsymbol{y}}_{n}=[\hat{\boldsymbol{y}}^{(1)}_{n},\dots,\hat{\boldsymbol{y}}^{(T)}_{n}]\in\hat{Y}$
where $\hat{Y}$ represents the set of all the outputs. We denote the label of
the $n$-th input data as
$\boldsymbol{y}_{n}=[\boldsymbol{y}^{(1)}_{n},\dots,\boldsymbol{y}^{(T)}_{n}]\in
Y$ where $Y$ represents the set of all the labels for the input data. The loss
function, denoted by $L$, is defined as the normalized squared difference
between labels and RNN’s output,

$L\left(\hat{Y},Y\right)=\frac{1}{|Y|}\sum_{\boldsymbol{y}_{n}\in
Y}(\hat{\boldsymbol{y}}_{n}-\boldsymbol{y}_{n})^{2}$ (6)

where the cardinality $|Y|$ is defined as the number of data instances in the
set $Y$. We use BTT to compute the derivatives of Eqn. (6) with respect to the
weights and bias vectors. We update the weights using the gradient descent
based method, called Adaptive Moment Estimation (ADAM) [24]. ADAM is an
iterative optimization algorithm that is widely used in modern machine
learning algorithms to minimize loss functions where it employs the averages
of both the first moment gradients and the second moment of the gradients for
computations. ADAM generally converges faster than standard gradient descent
methods and saves memory by not accumulating the intermediate weights.

##  III Performance analysis

In this section, we work on three diverse problems: 1) a Lorenz system having
a formulation error, 2) denoising of noisy trajectories of particle swarms in
collective motion, and 3) forecasting of groundwater streamflow of hydrologic
catchments. These three problems, in that order, represent three fields,
ordinary differential equations, collective motion, and hydrological modeling.

Here, we use the RNN model available in the machine learning library _Pytorch_
[25], which requires the parameter inputs sequential length, input size,
learning rate, output size, number of RNNs, and number of hidden nodes. The
sequential length is the length of the input data sequence, which is $T$ in
Fig. 2. The input size is the number of features or dimensions in the dataset,
which is the length of the vector $\boldsymbol{x}^{(t)}$ in Fig. 2. The
learning rate is associated with the ADAM optimization routing where a large
learning rate allows the RNN to learn faster at a cost of arriving into a sub-
optimal final set of weights. A smaller learning rate may allow the RNN to
learn further optimal or even globally optimal set of weights but may take
significantly longer time to train. The number of RNNs is the number of basic
RNN units that are stacked together which is equal to the number of hidden
layers, where it is $K=3$ in Fig. 2. The number of hidden nodes is the user
input number of nodes at each hidden layer of the RNN, which is $N_{h_{k}}$
for $k$-th hidden layer of Fig. 2. Finding both the optimal number of RNNs and
the optimal number of hidden nodes could prevent possible over-fitting or
under-fitting scenarios of the RNN’s training process.

###  III.1 Lorenz system with a formulation error

The Lorenz system is a chaotic dynamical system [26], that is formulated as a
set of ordinary differential equations which is known for having chaotic
solutions for certain parameter values and initial conditions [27]. Here, we
assume that we are given a Lorenz system with a formulation error so that our
task is to make an ANN to generate the true solutions related to the input
erroneous solutions generated by the system with the formulation error. The
standard Lorenz system of the temporal variable $t$ and three spacial
variables $y_{1}$, $y_{2}$, and $y_{3}$ is given by

$\displaystyle\frac{dy_{1}}{dt}=\sigma(y_{2}-y_{1}),$ (7)
$\displaystyle\frac{dy_{2}}{dt}=y_{1}(\rho-x_{3})-y_{2},$
$\displaystyle\frac{dy_{3}}{dt}=y_{1}y_{2}-\beta y_{3},$

where $\sigma$, $\rho$, and $\beta\in\mathbb{R}$ are parameters. To mimic a
formulation error, we eliminate the term $-y_{2}$ from Eqn. (7) as

$\displaystyle\frac{dx_{1}}{dt}=\sigma(x_{2}-x_{1}),$ (8)
$\displaystyle\frac{dx_{2}}{dt}=x_{1}(\rho-x_{3}),$
$\displaystyle\frac{dx_{3}}{dt}=x_{1}x_{2}-\beta x_{3}.$

Note that, we replaced the variable $y$ in Eqn. (7) with variable $x$ to
distinguish data from two systems. We set $\sigma=3$, $\rho=28$, and
$\beta=8/3$ in both Eqn. (7) and Eqn. (8) since the solutions with these
parameters are well-known for showing chaotic behaviors. For different initial
conditions, we solve the system in Eqn. (7) and treat them as training labels
for the RNN. Now, we generate the erroneous solutions for Eqn. (8) based on
Eqn. (7) to use as training data as we will explain in the sequel. We use the
trained RNN to generate the true solutions for the input erroneous solutions.

We generate 500 random initial conditions, $\boldsymbol{y}_{n}^{(0)}=$
$\left(y_{n,1}^{(0)},y_{n,2}^{(0)},y_{n,3}^{(0)}\right)$ where 1 $\leq$ $n$
$\leq$ 500, from the uniform distribution $\mathbb{U}[-15,15]$. We solve Eqn.
(7) and generate 500 orbits, denoted as $N$ with the step-size, denoted as
$\delta$, is 0.01, and the total time-steps, denoted as $T$, is 5000 (i.e.,
$t\in[0,50]$ of Eqn. (7)). Here, the $n$-th orbit is given as
$\boldsymbol{y}_{n}=\left[\boldsymbol{y}_{n}^{(0)}\right.$, $\dots$,
$\boldsymbol{y}_{n}^{(t)}$, $\dots$,
$\left.\boldsymbol{y}_{n}^{(5000)}\right]$, where
$\boldsymbol{y}_{n}^{(t)}={\left(y_{n,1}^{(t)},y_{n,2}^{(t)},y_{n,3}^{(t)}\right)}^{T}$.
Now, we explain the process of generating the erroneous orbit
$\boldsymbol{x}_{n}=\left[\boldsymbol{x}_{n}^{(0)}\right.$, $\dots$,
$\boldsymbol{x}_{n}^{(t)}$, $\dots$,
$\left.\boldsymbol{x}_{n}^{(5000)}\right]$ from the above true orbit. For some
$t$, we treat $\boldsymbol{y}_{n}^{(t)}$ in the $n$-th orbit as an initial
condition for the system in Eqn. (8) and solve it for some time-steps, denoted
as $\eta$, where $\eta<T$, with the same values for the parameters $\sigma$,
$\rho$, $\beta$, and $\delta$ as above. The parameter $\eta$ governs the
corruption level of the orbits such that a big value of $\eta$ contributes a
big level of corruption to the orbits. We treat the last point of the solution
set, that is the $\eta$-th point, as the erroneous solution relevant to the
correct point $\boldsymbol{y}_{n}^{(t)}$, that we denote by
$\boldsymbol{x}_{n}^{(t)}$, see Fig. 3. We carry out this process for each $n$
and generate the entire erroneous orbit
$\left[\boldsymbol{x}_{n}^{(0)}\right.$, $\dots$, $\boldsymbol{x}_{n}^{(t)}$,
$\dots$, $\left.\boldsymbol{x}_{n}^{(5000)}\right]$. Similarly, we compute the
erroneous orbits corresponding to each of the 500 true orbits that we
generated by changing initial conditions. The first 400 orbits of the system
in Eqn. (8) are chosen as the training set, and the corresponding 400 orbits
that are derived from the system in Eqn. (7) are chosen as their corresponding
labels.

Figure 3: For some $t^{\prime}$, where $1\leq t^{\prime}\leq T$, we treat
$\boldsymbol{y}_{n}^{(t^{\prime})}$ in the $n$-th true orbit as the initial
condition for the erroneous Lorenz system and solve it for some time-steps,
denoted as $\eta$, with the same parameter values used for solving the true
Lorenz system. We consider the last point $\eta$, denoted by
$\boldsymbol{x}_{n}^{(t^{\prime})}$, of the solution set as the erroneous
solution relevant to the true point $\boldsymbol{y}_{n}^{(t^{\prime})}$. We
carry out this process for each $t^{\prime}$ to generate the corrupted orbit
$\left[\boldsymbol{x}_{n}^{(0)}\right.$, $\dots$, $\boldsymbol{x}_{n}^{(t)}$,
$\dots$, $\left.\boldsymbol{x}_{n}^{(T)}\right]$. Here, $\eta$ can be
considered as the corruption level of the orbit where a big $\eta$ contributes
a high corruption.

First, we set $\eta=1$ and generate 500 true and erroneous orbits each with
the 5000 time-steps as explained above, where the first 400 orbits are used
for training while the last 100 orbits are used for testing. In order to avoid
both the over-fitting and the semi-convergence of RNN, the best number of
training epochs was chosen by running the RNN twice as follows: first, we run
the RNN with a fixed, but big, number of epochs to find the epoch number that
gives the minimum loss; second, we re-run the RNN for that many epochs to
complete the training. For the Lorenz system, based on the corruption level
$\eta$, the optimum number of epochs ranges between 179 and 930. We utilize
the activation function ReLU as it improves the reconstruction error about 10
times than that of the hyperbolic tangent activation function. In order to
find the optimal parameters, which are learning rate, number of RNN, number of
hidden node; first, we discretize the parameter space; then, we train and test
the RNN with those discretized parameter combinations; finally, we select the
best values for the parameters that give the least reconstruction error. In
all the experiments, a data batch is equal to the entire dataset. We set the
other parameters of RNN as follows: sequential length to 5000 since each orbit
has 5000 time-steps; input size to three since each orbit is three-
dimensional; output size to three since the output is also three-dimensional;
the learning rate of the Lorenz orbits to 0.01.

We set the RNN to three hidden recurrent layers with $256$ nodes each. For
each $1\leq n\leq 400$, we input each $\boldsymbol{x}_{n}$ into the RNN one at
a time and generate the corresponding output of the RNN, denoted as
$\hat{y}_{n}$. We compute the reconstruction error $L$ between the outputs and
the labels using Eqn. (6). We use BTT to compute the derivatives of Eqn. (6)
with respect to the weights and the bias vectors, and update the weights and
the bias vectors using ADAM as explained in Sec. II.2. The remaining 100 test
orbits, $\boldsymbol{x}_{n}$ where $401\leq n\leq 500$, are fed into this
trained RNN one orbit at a time and obtain the recovered solutions,
$\hat{\boldsymbol{y}}_{n}^{(t)};401\leq n\leq 500$. We repeat the same
experiment nine more times with $\eta=2,3\dots,9$. Fig. 4 shows two-
dimensional projections of the three-dimensional true, erroneous, and
recovered orbits for the cases $\eta$ = 1, 5, and 9.

Figure 4: Recovery of the true solution of an erroneous Lorenz system using
RNNs. Here, all the figures are two-dimensional ($y$-$z$ plan) projections of
the three-dimensional orbits. 120 true orbits are generated using a Lorenz
system and a Lorenz system with a formulation error is used to generate
corrupted orbits corresponding the true orbits with some known corruption
level $\eta$. Here, while (a-c) show three arbitrarily chosen true orbits,
(d-f) show the erroneous orbits corresponding to them, respectively, where the
green disks represent the starting point and red disks represent the ends. For
each $\eta$, we set the parameters of the RNN as follows: 5000 for sequential
length, 3 for both input size and output size, 0.01 for learning rate, three-
stacked RNN with 128 number of hidden node each, and ReLu for the activation
function. Moreover, the best number of epochs vary between 179 and 930 based
on $\eta$. First, we train an RNN with arbitrary 100 orbits for each $\eta$,
and then we pass the erroneous orbits through the RNNs and recovery them as
show in (g)-(i). The start and the end points of the recovered orbits are
close to those of the true orbits whereas they are far apart between erroneous
and true orbits. Root mean square error of the reconstructions are 0.072,
0.15, 0.27 for the corruption levels $\eta=1,5$, and 9, respectively.

We observe that due to the formulation error of the Lorenz system the
erroneous orbits show some shifting of the initial and the end points from
their true locations. This shift does not only occur at the two extremes but
all the intermediate points are also shifted along the orbit in the same
direction. We observed that the recovered orbits have their extreme points
around their true locations. Moreover, the recovery quality gets worse as
$\eta$ increases. To analyze the recovery performance with increasing
corruption levels, we compute the normalized root mean square error, denoted
by RMSE, between the recovered and the true orbits as

$RMSE=\sqrt{\frac{1}{|Y|}\sum_{\forall\boldsymbol{y}_{n}\in
Y}\|\boldsymbol{y}_{n}-\hat{\boldsymbol{y}}_{n}\|_{F}^{2}},$ (9)

where $|Y|$ denotes the carnality, i.e, the number of observations in the test
set, and $F$ denotes the Frobenius norm. Fig. 5 shows that the RMSE of the
recovery with respect to the corruption level where we see that the
reconstruction error increases rapidly with respect to the increase of the
corruption level.

Figure 5: RMSE of the reconstruction of corrupted Lorenz orbits versus
different levels of the corruptions ($\eta$). 120 true orbits are constructed
using a conventional Lorenz system by changing initial conditions. Then, these
true orbits are used as the initial conditions to generate solutions for the
Lorenz system with a formulation error where this generates 120 corrupted
orbits. Each of these corrupted orbits has undergone $\eta$-level of
corruption where $\eta=1,2,\dots,9$ to create 9 separate datasets each with
120 true orbits and 120 corrupted orbits. Then, for each $\eta$, an RNN having
three hidden layers is trained with the first 100 orbits and then test on the
remaining 20 orbits. The root mean squared errors (RMSE) between the true and
the reconstructed trajectories justify that the reconstruction error increases
rapidly when the corruption level increases.

###  III.2 Noisy trajectories of collective motion

Analysis of the trajectories of collectively moving agents such as fish [28]
and birds [29], is a highly active field in computer vision. Tracked
trajectories of such agents are often contaminated with noise due to the
causes like insufficient camera precision and lack of accuracy of the tracking
method [30]. Here, we use RNNs to reconstruct the true collective motion
trajectories from observed noisy trajectories. We simulate the collective
motion trajectories using a modified version of a classic self-propelled
particle model, named the Vicsec model [13].

The collective motion is defined as a spontaneous emergence of the ordered
movement in a system consisting of many self-propelled agents. The Vicsek
model given in Ref. 13 is one of the widely used models that describe such
behavior. This model performs simulations on a square-shaped region with
periodic boundary conditions. We generate a synthetic collective motion
dataset by using the generalized version presented in Ref. 12 of the classic
Vicsek model where we incorporated a rotational matrix, denoted as
$R^{(t)}_{n}$. The agent-wise temporal rotation matrix $R^{(t)}_{n}$ imposed
to the $n$-th agent at the $t$-th time-step allows us to simulate interesting
collective motion scenarios while ensuring the intra-group interactions
between the agents. Based on the average direction of motion of all the
particles in a neighborhood, denoted by $N^{(t)}_{n}$, within a radius $r_{d}$
of the $n$-th particle, the generalized Vicsek model updates the orientation
of the $n$-th particle at the $t$-th time-step, denoted by
$\theta_{n}^{(t)}\in[-\pi,\pi)$. Therefore, the orientation of the $n$-th
agent at the $t$-th time-step, defined as $\theta_{n}^{(t+1)}$, is computed as

$\theta_{n}^{(t+1)}=\arg(\boldsymbol{u}^{(t)}_{n})+\epsilon_{n}^{(t)},$ (10)

where, $\epsilon_{n}^{(t)}$ is a noise parameter imposed to the orientation of
the $n$-th agent at the $t$-th time-step, and

$\boldsymbol{u}^{(t)}_{n}=\frac{1}{|N^{(t)}_{n}|}\sum_{j\in
N^{(t)}_{n}}R^{(t)}_{j}\begin{bmatrix}\cos(\theta_{j}^{(t)})\\\
\sin(\theta_{j}^{(t)})\end{bmatrix}.$ (11)

The position vector of the $n$-th agent at the $t$-th time-step is given as

$\boldsymbol{x}_{n}^{(t+1)}=\boldsymbol{x}_{n}^{(t)}+v_{n}^{(t)}R_{n}^{(t)}\begin{bmatrix}\cos(\theta_{n}^{(t)})\\\
\sin(\theta_{n}^{(t)})\end{bmatrix}\delta,\\\ $ (12)

where, $v_{n}^{(t)}$ denotes the speed of the $n$-th agent at the $t$-th time-
step and $\delta$ is the step size. Eqns. (10), (11), and (12) with
$R_{n}^{(t)}=[1,0;0,1]$ for all $n$’s and $t$’s implies the classic Vicsek
model.

We use this rotational matrix to formulate a spiral collective motion scenario
since spiral collective motion is often observed in nature ranging from
biology, such as bacteria colonies [31], to astronomy, such as spiral galaxies
[32]. Here, first, we formulate an anticlockwise Archimedean spiral that
rotates an angle of $3\pi$ and then compute the agent-wise temporal rotational
matrix based on that spiral. We assume that the two-dimensional coordinates
$\left[\boldsymbol{c}^{(1)};\dots;\boldsymbol{c}^{(t)};\dots;\boldsymbol{c}^{(T)}\right]$
where $\boldsymbol{c}^{(t)}\in\mathbb{R}^{2}$ represent this spiral such that

$\begin{split}\boldsymbol{c}^{(t)}=r^{(t)}\begin{pmatrix}\cos\left(\kappa^{(t)}\right)\\\
\sin\left(\kappa^{(t)}\right)\end{pmatrix}\ ;\ \text{where}\\\
r^{(t)}=1+\frac{3}{T-1}(t-1)\ \text{and}\
\kappa^{(t)}=\frac{3\pi}{T}(t-1),\end{split}$ (13)

for $t=1,\dots,T$. Here, $r^{(t)}$ is the variable radius of the spiral where
it changes from 1 to 4 and $\kappa^{(t)}$ is the angle of rotation with
respect to the origin that varies from 0 to $3\pi$. The rotation angle with
respect to the $(t-1)$-th coordinates,
$\boldsymbol{c}^{(t-1)}=\left[c^{(t-1)}_{1},c^{(t-1)}_{2}\right]^{T_{r}}$, of
this spiral, denoted by $\gamma^{(t)}$, is

$\gamma^{(t)}=\tan^{-1}\left(\frac{c_{2}^{(t)}-c_{2}^{(t-1)}}{c_{1}^{(t)}-c_{1}^{(t-1)}}\right)\
;\ t=2,\dots,T.$ (14)

Thus, the two-dimensional rotational matrix for the $i$-th agent at the $t$-th
time-step is given by

$R^{(t)}_{n}=\begin{pmatrix}\cos\left(-\gamma^{(t)}\right)&-\sin\left(-\gamma^{(t)}\right)\\\
\sin\left(-\gamma^{(t)}\right)&\cos\left(-\gamma^{(t)}\right)\end{pmatrix}$
(15)

Eqns. (10), (11), and (12) with the rotational matrix in Eqn. (15) provide the
complete formulation of the system generating the spiral collective motion
dataset.

We generate 30 agents using Eqns. (10), (11), (12), and (15) with 201 for the
time-steps ($T$), a rectangular domain with periodic boundary conditions. We
set two for the radius of interaction ($r_{d}$), 0.05 for the speed of the
particles ($\boldsymbol{v}^{(t)}_{i}$ for all $t$ and $i$), 0.05 for the noise
on the orientation ($\epsilon^{(t)}_{i}$ for all $t$ and $i$), and one for the
time-step size ($\delta$), see the first row of Fig. 6. Then, we impose three
noise levels sampled from the Gaussian distribution
$\mathcal{N}(0,\sigma^{2})$, where $\sigma$ = .2, .4, and .6, into the
original dataset and make three copies of that, see the second row of Fig. 6.

Similar to the Lorenz experiment, in order to avoid both the over-fitting and
the semi-convergence of RNNs, the best number of training epochs was chosen by
training an RNN with a fixed, but big, number of epochs to find the epoch
number that gives the minimum loss, and then training another RNN with that
many epochs. Based on the noise level of the collective motion trajectories,
the best number of epochs is between 6548 and 14032. We tested both the
hyperbolic tangent and Relu as activation functions and found that the loss
with hyperbolic tangent is smaller than the one with ReLu. Thus, we set the
activation function of the RNN to a hyperbolic tangent. Similar to Lorenz
experiment, to find the optimal parameter values for the learning rate, the
number of RNNs, the number of hidden nodes, first, we discretize the parameter
space, then, we train and test the RNN with those discretized parameter
combinations, finally, we select the best values for the parameters that give
the least reconstruction error. Thus, the best learning rate is 0.0005, the
best number of RNNs is two, and the best number of nodes is 64. Moreover, the
other parameters of the RNN are set as follows: a batch to one trajectory,
sequential length to 201 since each trajectory has 201 time-steps, input size
to two since each trajectory is two-dimensional, output size to two since the
output is also two-dimensional.

For each noise level, we train one RNN with 80%, i.e., 24, of the
trajectories, and test with the rest of them. We input each trajectory
$\boldsymbol{x}_{n}$ into the RNN one at a time and generate the corresponding
output of RNN, denoted as $\hat{\boldsymbol{y}}_{n}$. We compute the
reconstruction error between $\hat{\boldsymbol{y}}_{n}$ and the label
$\boldsymbol{y}_{n}$ using Eqn. (9). We use BTT to compute the derivatives of
Eqn. (6) with respect to the weights and update the weights using ADAM as
explained in Sec. II.2. The remaining 20%, i.e., six, of the noisy
trajectories are fed into the corresponding trained RNN one trajectory at a
time and obtain the reconstructions, $\hat{\boldsymbol{y}}_{n}$, see the third
row of Fig. 6. Root mean square errors of the reconstructions, computed using
Eqn. (9), are $0.05,0.08,0.10$ for the noise levels $\sigma$ = .2, .4, and .6
of $\mathcal{N}(0,\sigma^{2})$, respectively. This test justifies that the
RNN’s performance decreases as the noise level increases.

Figure 6: Reconstruction of noisy collective motion trajectories using RNNs.
(a-c) An original trajectory dataset of 30 agents is generated using a
generalized version of the Vicsek model (all the three datasets in a-c are
identical). We make three copies of the original dataset by corrupting it with
three levels of Gaussian noises with $\sigma$ = .2, .4, and .6 of
$\mathcal{N}(0,\sigma^{2})$ where (d-f) show 20% of the noisy trajectories
chosen for the testing. 80% of the trajectories of each dataset are used for
the training of the RNNs, while the rest is used for the testing. We set the
parameters of the RNNs as follows: 201 for the sequential length, two for both
the input size and the output size, 0.0005 for the learning rate, two-stacked
RNN with 64 number of hidden nodes in each, 6548-14032 for the number of
epochs based on the noise level, tanh for the activation function. We train
each RNN with the noisy trajectories as the inputs and the original
trajectories as the outputs. Then, we pass the noisy test trajectories in each
(d-f) through the corresponding RNN and denoise them as shown in (g)-(i) where
the denoised trajectories are shown in continuous lines while the
corresponding original trajectories are shown in dashed lines with the same
colors. Root mean square errors of the reconstructions are $0.05,0.08,0.10$
for the noise levels $\sigma$ = .2, .4, and .6 of $\mathcal{N}(0,\sigma^{2})$,
respectively.

###  III.3 Spiky time series of rainfall-runoff

We use RNNs to forecast on hydrology data, especially rainfall-runoff time
series which tends to be spiky. Then, compare the RNN’s results with another
forecast produced by a widely used hydrological model named as GR4J available
in Ref. 33. Hydrology data, especially, rainfall-runoff, is highly volatile
with frequent spikes so that the hydrologists have been struggling over the
past decades to improve modeling of them [34]. GR4J is a daily lumped
rainfall-runoff model proposed to understand hydrologic catchments’ behaviors
[33]. This model is characterized by two independent variables, namely,
precipitation (denoted by $p$) and potential evapotranspiration (denoted by
PET), dependent variable, namely, streamflow (denoted by $q$), and four
parameters, namely, the capacity of the production store (denoted by
$\nu_{1}$), ground exchange coefficient (denoted by $\nu_{2}$), the capacity
of the nonlinear routing store (denoted by $\nu_{3}$), and unit hydrograph
time base (denoted by $\nu_{4}$). The GR4J model was used in combination with
the degree-day snowmelt module available in Ref. 35, to account for snow in
northern latitudes. This degree-day snowmelt module has four parameters,
namely, threshold temperature (denoted by TT), degree-day factor (denoted by
CFMAX), refreezing factor (denoted by CFR), and water holding capacity of the
snowpack (denoted by CWH). GR4J requires tuned values for eight parameters in
order to produce accurate forecasts. We use 10 years of recent precipitation,
evapotranspiration, and streamflow data from three arbitrarily chosen U.S.
catchments, that we denote as Site 1 (MOPEX ID is 2365500), Site 2 (MOPEX ID
is 2492000), and Site 3 (MOPEX ID is 7378000), of the model parameter
estimation Experiment (MOPEX) database [36].

From the 10 recent years, we use the first seven years, i.e., 2557 days, for
the training and the last three years, i.e., 1096 days, for the forecasting.
For all three hydrological sites, we follow the same training and forecasting
routings and use the same parameter values in RNNs. Similar to the previous
experiments, we utilize hyperbolic tangent as the activation function and
0.001 as the learning rate. We use three-stacked RNNs with 512 nodes at each
hidden layer. We utilize a special training and forecasting approach here
using sliding windows as it aids sampling of adequate amount of training data
just from the given time series. We make overlapping windows of the data such
that each window is $L$ days long for some $L<2557$ where two consecutive
windows are misaligned by only one day, see Fig. 7 for a special case of
$L=5$. Since we vary window size in this experiment such that
$L=[5,15,30,45,60,10,240,365]$, we train one RNN for each case of $L$. Thus,
the sequential length parameter of the RNNs is set to the corresponding value
for $L$ one at a time. Since RNN’s input is two one-dimensional input vectors
and the output is a single one-dimensional vector, the input size parameter is
set to two and the output size parameter is set to one. Similar to the
previous experiments, in order to avoid both the over-fitting and the semi-
convergence of the RNN, the best number of training epochs was chosen by
running the RNN with a fixed, but big, number of epochs to find the epoch
number that gives the minimum loss, and then re-running the RNN with that many
epochs to get the results. Based on $L$, the best number of epochs is between
99178 and 99994. We set 0.001 for the learning rate, 99200 for the number of
epochs, tanh for the activation function.

We denote $\boldsymbol{x}_{t}=[p,PET,q]$ for $t\in[1,3653]$ where $p$ and PET
are used as the inputs to the RNN and $q$ is used as the labels during the
training process. Specifically, for the first iteration, we train the RNN with
$X_{1}=[p;PET]_{2\times T}$ as the input and $Y_{1}=[q]_{1\times T}$ as the
labels. Similarly, for all $t\in[2,2558-L]$, we train the same RNN with all
such $X_{t}$’s and $Y_{t}$’s, see the blue color pairs in Fig. 7. For
$t\in[1,3654-L]$, we pass each $X_{t}$ through the RNN and obtain the modeled
streamflow windows $Y_{t}$, see the blue-green color pairs in Fig. 7. For each
$t\in[1,3653]$, the corresponding streamflow is computed as the average of the
streamflows between overlapped windows. The modeled streamflow for
$t\in[1,2557]$ is used to validate the training performance of the RNN and the
modeled streamflow for $t\in[2558,3653]$ is considered as the streamflow
forecast, see the red color vectors in Fig. 7.

Figure 7: Training and forecasting RNNs with recent 10 years of MOPEX data
using a sliding window approach. We train an RNN with sliding windows of
precipitation ($p$), potential evapotranspiration (PET), and streamflow ($q$)
where the first seven years, i.e., 2557 days, are used for the training to
make forecasting for the last three years, i.e., 1096. We make overlapping
windows of the data such that each window is $L$ time-steps long for some
$L<2557$ where this figure shows a special case of $L=5$. Here, each
$\boldsymbol{x}_{t}$ for $t\in[1,3653]$, is a vector of $p$, PET, and $q$ so
that $p$ and PET are used as the inputs to the RNN and $q$ is used as the
labels during the training process, see the blue color vector pairs.
Specifically, for the first iteration, we train the RNN with
$X_{1}=[p;PET]_{2\times T}$ as the input and $Y_{1}=[q]_{1\times T}$ as the
labels. Similarly, for all $t\in[2,2558-L]$, we train the same RNN with all
such $X_{t}$’s and $Y_{t}$’s. For $t\in[1,3654-L]$, we pass each $X_{t}$
through the RNN and obtain the modeled streamflow windows $Y_{t}$, see the
blue-green color pairs. For each $t\in[1,3653]$, the corresponding streamflow
is computed as the average of the streamflows between overlapped windows. The
modeled streamflow for $t\in[1,2557]$ is used to validate the training
performance of the RNN and the modeled streamflow for $t\in[2558,3653]$ is
considered as the streamflow forecast, see the red color.

We compute RMSEs between the true and the forecasted streamflows using Eqn.
(9) for each widow size $L$. Fig. 8 shows RMSEs versus window size for Site 1
where we observe the least RMSE at $L=45$. The orange color plots of Figs.
9(a), (c), and (e) show the true streamflows and the red color plots of those
figures show the modeled streamflow by RNNs, and the red color plots in Figs.
9(b), (d), and (f) show the absolute difference between the true and the RNN’s
modeled streamflows. The RMSEs, computed using Eqn. (9), between the true and
forecasted streamflows are 0.9, 1.1, and 2.0 for Sites 1, 2, and 3,
respectively.

Figure 8: RMSE of Site 1 (MOPEX ID is 2365500) for varying window sizes ($L$)
where the RMSE is computed between the true streamflow and RNN’s forecasted
streamflow. We observe that the least RMSE is attained around $L=50$.

The eight parameters of GR4J have the recommended ranges as follows:
$\nu_{1}\in[0,1000]$, $\nu_{2}\in[-5,5]$, $\nu_{3}\in[0,300]$,
$\nu_{4}\in[0.5,5]$, TT $\in[-3,3]$ , CFMAX $\in[0,20]$, CFR $\in[0,1]$, and
CWH $\in[0.0.8]$ [35, 33, 34]. We uniformly discretize GR4J’s parameter space
of the eight parameters into 50000 points where each point represents a
combination of parameter values. We run GR4J with all of the parameter
combinations for first seven years of the data where the best parameter
combination with respect to RMSE is $\nu_{1}$ = 571, $\nu_{2}$ = -0.03,
$\nu_{3}$ = 48, $\nu_{4}$ = 2, TT = -0.20, CFMAX = 4.41, CFR = 0.36, and CWH =
0.23, for the Site 1. We set the same values for the parameters of Site 2 and
Site 3 since their optimum parameter values are either the same or
sufficiently close to those of Site 1. We calibrate the GR4J model with the
data for each site and model the streamflow $q_{t}$ for $t\in[1,3653]$ that we
show by blue color plots in Figs. 9 (a), (b), and (c). The absolute difference
of the true and GR4J’s modeled streamflows are shown by the blue color plots
in Figs. 9 (b), (d), and (f). The RMSEs between the true and the forecasted
streamflows are 1.0, 1.3, and 2.7 for the sites 1, 2, and 3, respectively. We
observe that the RMSEs associated with the RNN’s modeled streamflows are less
than that of the Gr4J’s modeled streamflows. GR4J is a widely used streamflow
modeling scheme[37]; however, the performance of RNN in forecasting of
streamflow is capable of exceeding the performance of GR4J.

Figure 9: Comparison of streamflow, denoted by $q=[q_{t}|\ \forall t]$,
forecasting performance between RNN and GR4J for three hydrological sites,
named as Site 1 (MOPEX ID is 2365500), Site 2 (MOPEX ID is 2492000), and Site
3 (MOPEX ID is 7378000). We use 10 recent years of data for the study where
the first seven years, i.e., 2557 days, are used for training and the last
three years, i.e., 1096 days, are used for the testing. We set the same value
for the same parameter between each of the RNN implemented for each site as
follows: 45 for the sequential length, 2 for the input size, 1 for the output
size, 0.001 for the learning rate, three-stacked RNNs, 512 for the number of
hidden nodes, 99200 for the number of epochs, tanh for the activation
function. We train the RNN with the first seven years of data where the inputs
are precipitation, denoted by $p$, and potential evapotranspiration, denoted
by PET, and the output is streamflow. Then, we pass [$p$; PET] through the
trained RNN to obtain the modeled $q$ where $q$ for $t\in[1,2557]$ is used for
training validation and $q$ for $t\in[2558,3653]$ is considered as the
forecasting. Orange color plots in (a), (c), and (e) show true streamflow, and
the red color plots therein show the RNN’s modeled streamflow where the
vertical black color line separates training and testing periods. Since the
red color is not visible since it is underneath the orange color, we compute
the absolute difference between the true and the modeled streamflows that we
show by red color plots in (b), (d), and (f). We set the parameters of GR4J,
namely, capacity of the production store ($\nu_{1}$), ground exchange
coefficient ($\nu_{2}$), the capacity of the nonlinear routing store
($\nu_{3}$), and unit hydrograph time base ($\nu_{4}$), threshold temperature
(TT), degree-day factor (CFMAX), refreezing factor (CFR), and water holding
capacity of the snowpack (CWH), to 571, -0.03, 48, 2, -0.02, 4.41, 0.36, and
0.23, respectively, for all the three sites and calibrate it. Then, we model
the streamflow similar to that of RNN for the entire period where the blue
color plots in (a), (c), and (e) show the GR4J’s modeled streamflows and the
blue color plots in (b), (d), and (f) show the absolute differences between
true and modeled streamflows. The RMSE between forecasted and observed
streamflows of the Sites 1, 2, and 3 are 0.9, 1.1, and 2.0 for RNNs,
respectively, and are 1.0, 1.3, and 2.7 for GR4J, respectively.

Table 2 summarizes the parameter values used in RNNs for each of the dataset.

Table 2: Parameter values used in the RNNs trained for each datatset. Dataset |  Sequantial length |  Input size |  Output size |  Learning rate |  Num. of RNNs |  Num. of hidden nodes |  Num. of epochs |  Activation function
---|---|---|---|---|---|---|---|---
Lorenz system with a formulation error |  5000 |  3 |  3 |  0.01 |  3 |  128 |  179-930 |  ReLu
Noised trajectories of collective motion |  201 |  2 |  2 |  0.0005 |  2 |  64 |  6548-14032 |  tanh
Spiky time series of rainfall-runoff |  5–365 |  2 |  1 |  0.001 |  4 |  512 |  99178-99994 |  tanh

##  IV Conclusion

In this paper, we have utilized RNNs for three diverse tasks, namely,
correction of a formulation error, reconstruction of corrupted particle
trajectories, and forecasting of streamflow data, in three diverse fields of
dynamical systems, namely, ODEs, collective motion, and hydrological modeling,
respectively. The traditional approaches to solve spatiotemporal dynamical
systems may not capture non-linear and complex relationships in the data since
they are mostly either linear or non-linear model-based. Such spatiotemporal
nonlinear systems can effectively be learned by nonlinear and model-free
machine learning models such as RNNs with both high robustness and better
accuracy as RNNs possess internal memory that enhance the learning process.

We used the ODE model Lorenz system to produce two sets of data such that one
consist the solutions for a correct model whereas the other consist the
solution for an erroneous Lorenz model. We observed that the trained RNN was
capable of producing the corresponding correct solution when the incorrect
solution was provided. This approach can be used for eliminating errors
associated with ODEs solvers even when only an imprecise closed-form solution
is available. Future work in this context maybe using a further advanced
neural network tool such as LSTM presented in Ref. 8 as it is composed of
three gates that regulate a better information flow through the unit. Such an
approach enhances the learning process of spatiotemporal data so that it
better aid in eliminating the formulation error even for longer corrupted
orbits.

We used a generalized version of the classic Vicsek model to generate
trajectories imitating a spiral collective motion scenario. We diagnosed the
influence of noise contamination on RNN’s denoising performance. Such noise
contamination in collective motion trajectories can occur due to the causes
like less precision in video cameras and less robustness of the tracking
methods [38]. The results show that RNNs are capable of denoising the noisy
trajectories while ensuring the pattern of the underlying collective motion.
Such an ANN-based denoising method can be integrated into multi-object
tracking methods to generate trustworthy results. As future work, we are
planning to use RNN’s denoised trajectories for the subsequent collective
motion analysis tasks such as phase transition detection as we presented in
Ref. 39.

Streamflow modeling is a sophisticated task in the field of hydrology due to
its high volatility and frequent spikes. GR4J is a popular tool for streamflow
modeling since it has shown trustworthy streamflow modeling performance.
However, we use RNNs to model streamflow with a sliding window approach. The
sliding window approach is capable of generating a big amount of homogeneous
data from a given time series. Moreover, the sliding window approach can be
used to find the most effective sequential length of the data samples that
captures the most features of the given time series. Thus, we trained an RNN
with sliding windows of the empirically tested best length and then model
streamflow for the entire time duration. We observed that the RNN’s training
process involves less reconstruction error than that of GR4J. The reason for
that is RNN is model-free whereas GR4J is model-based, so RNN has more
tendency to adjust to any natural system than that of Gr4J. The observation
that the streamflow forecasting of RNN is better than that of GR4J justifies
the best performance of RNN in real-life time series forecasting. We are
planning to model streamflow using other well-known models such as HBV [40]
and Simhyd [41], and then compare those results with that of RNNs.

RNNs are capable of learning spatiotemporal characteristics of data derived
from dynamical systems that we used for three such applications representing
three fields. The observations based on our broad analysis justify that an RNN
is an effective machine learning tool in learning a vast array of dynamical
systems. The internal memory module in RNNs grants this ability which is not
available in most of the spatiotemporal dynamic analyzing tools.

## Acknowledgments

The authors would like to thank the Google Cloud Platform for granting
Research Credit to access its GPU computing resources under the project number
397744870419. E. M. Bollt has received funding from the Army Research Office
(No. N68164-EG), the National Institutes of Health (NIH-CRCNS) and also the
Defense Advanced Research Projects Agency.

## References

## References

  * Hu and Li [2016] M. Hu and F. Li, “A new method to solve numeric solution of nonlinear dynamic system,” Mathematical Problems in Engineering 2016 (2016).
  * Hornik, Stinchcombe, and White [1989] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural networks 2, 359–366 (1989).
  * Wang and Blum [1992] X. Wang and E. K. Blum, “Discrete-time versus continuous-time models of neural networks,” Journal of Computer and System sciences 45, 1–19 (1992).
  * Jayawardena [2013] A. W. Jayawardena,  _Environmental and hydrological systems modelling_ (CRC Press, 2013).
  * Diaconescu [2008] E. Diaconescu, “The use of NARX neural networks to predict chaotic time series,” Wseas Transactions on Computer Research 3, 182–191 (2008).
  * Yuan  _et al._ [2019] X. Yuan, P. He, Q. Zhu, and X. Li, “Adversarial examples: Attacks and defenses for deep learning,” IEEE Transactions on Neural Networks and Learning Systems 30, 2805–2824 (2019).
  * Sagheer and Kotb [2019] A. Sagheer and M. Kotb, “Unsupervised pre-training of a deep LSTM-based stacked autoencoder for multivariate time series forecasting problems,” Scientific Reports 9, 1–16 (2019).
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation 9, 1735–1780 (1997).
  * Cho  _et al._ [2014] K. Cho, B. Van Merriënboer, D. Bahdanau, and Y. Bengio, “On the properties of neural machine translation: Encoder-decoder approaches,” arXiv preprint arXiv:1409.1259 (2014).
  * Bianchi  _et al._ [2017] F. M. Bianchi, E. Maiorino, M. C. Kampffmeyer, A. Rizzi, and R. Jenssen, “An overview and comparative analysis of recurrent neural networks for short term load forecasting,” arXiv preprint arXiv:1705.04378 (2017).
  * Tabor [1989] M. Tabor,  _Chaos and integrability in nonlinear dynamics: an introduction_ (Wiley-Interscience, 1989).
  * Gajamannage and Paffenroth [2019] K. Gajamannage and R. Paffenroth, “Reconstruction of agents’ corrupted trajectories of collective motion using low-rank matrix completion,” in  _Proceedings of the 2019 IEEE International Conference on Big Data_ (IEEE, 2019) pp. 2826–2834.
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of self-driven particles,” Physical Review Letters 75, 1226–1229 (1995).
  * Csáji and Others [2001] B. C. Csáji and Others, “Approximation with artificial neural networks,” Faculty of Sciences, Etvs Lornd University, Hungary 24, 7 (2001).
  * Zhang  _et al._ [2016] S. Zhang, Y. Wu, T. Che, Z. Lin, R. Memisevic, R. R. Salakhutdinov, and Y. Bengio, “Architectural complexity measures of recurrent neural networks,” Advances in Neural Information Processing Systems 29, 1822–1830 (2016).
  * Bengio, Simard, and Frasconi [1994] Y. Bengio, P. Simard, and P. Frasconi, “Learning long-term dependencies with gradient descent is difficult,” IEEE Transactions on Neural Networks 5, 157–166 (1994).
  * Chollet [2017] F. Chollet,  _Deep Learning with Python_ (Manning Publications Co. LLC, chollet2017, 2017).
  * Siegelmann and Sontag [1991] H. T. Siegelmann and E. D. Sontag, “Turing computability with neural nets,” Applied Mathematics Letters 4, 77–80 (1991).
  * Schmidhuber [1992] J. Schmidhuber, “Learning to control fast-weight memories: An alternative to dynamic recurrent networks,” Neural Computation 4, 131–139 (1992).
  * Schäfer and Zimmermann [2007] A. M. Schäfer and H. G. Zimmermann, “Recurrent neural networks are universal approximators,” International Journal of Neural Systems 17, 253–263 (2007).
  * Werbos [1990] P. J. Werbos, “Backpropagation through time: what it does and how to do it,” Proceedings of the IEEE 78, 1550–1560 (1990).
  * Gruslys  _et al._ [2016] A. Gruslys, R. Munos, I. Danihelka, M. Lanctot, and A. Graves, “Memory-efficient backpropagation through time,” Advances in Neural Information Processing Systems 29, 4125–4133 (2016).
  * Manneschi and Vasilaki [2020] L. Manneschi and E. Vasilaki, “An alternative to backpropagation through time,” Nature Machine Intelligence 2, 155–156 (2020).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, “Adam: a method for stochastic optimization,” arXiv preprint arXiv:1412.6980 (2014).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” Advances in Neural Information Processing Systems 32, 8026–8037 (2019).
  * Devaney [2018] R. L. Devaney,  _A first course in chaotic dynamical systems: theory and experiment_ (CRC Press, 2018).
  * Lorenz [1963] E. N. Lorenz, “Deterministic nonperiodic flow,” Journal of Atmospheric Sciences 20, 130–141 (1963).
  * Gajamannage  _et al._ [2015a] K. Gajamannage, S. Butail, M. Porfiri, and E. M. Bollt, “Dimensionality reduction of collective motion by principal manifolds,” Physica D: Nonlinear Phenomena 291, 62–73 (2015a).
  * Gajamannage  _et al._ [2015b] K. Gajamannage, S. Butail, M. Porfiri, and E. M. Bollt, “Identifying manifolds underlying group motion in vicsek agents,” The European Physical Journal Special Topics 224, 3245–3256 (2015b).
  * Gajamannage  _et al._ [2021] K. Gajamannage, Y. Park, R. Paffenroth, and A. P. Jayasumana, “Reconstruction of fragmented trajectories of collective motion using hadamard deep autoencoders,” arXiv preprint arXiv:2110.10428 (2021).
  * Koizumi  _et al._ [2020] R. Koizumi, T. Turiv, M. M. Genkin, R. J. Lastowski, H. Yu, I. Chaganava, Q.-H. Wei, I. S. Aranson, and O. D. Lavrentovich, “Control of microswimmers by spiral nematic vortices: Transition from individual to collective motion and contraction, expansion, and stable circulation of bacterial swirls,” Physical Review Research 2, 033060 (2020).
  * Parnovsky and Parnowski [2010] S. L. Parnovsky and A. S. Parnowski, “Large-scale collective motion of RFGC galaxies,” Astrophysics and Space Science 325, 163–175 (2010).
  * Vaze  _et al._ [2011] J. Vaze, P. Jordan, R. Beecham, A. Frost, G. Summerell, J. Vaze, P. Jordan, R. Beecham, A. Frost, and G. Summerell, “Guidelines for rainfall-runoff modelling towards best practice model application. ewater cooperative research center,” Interim Version 1 (2011).
  * Jayathilake and Smith [2021] D. I. Jayathilake and T. Smith, “Assessing the impact of pet estimation methods on hydrologic model performance,” Hydrology Research 52, 373–388 (2021).
  * Kollat, Reed, and Wagener [2012] J. Kollat, P. Reed, and T. Wagener, “When are multiobjective calibration trade-offs in hydrologic models meaningful?” Water Resources Research 48 (2012).
  * Duan  _et al._ [2006] Q. Duan, J. Schaake, V. Andréassian, S. Franks, G. Goteti, H. Gupta, Y. Gusev, F. Habets, A. Hall, L. Hay, _et al._ , “Model parameter estimation experiment (mopex): An overview of science strategy and major results from the second and third workshops,” Journal of Hydrology 320, 3–17 (2006).
  * Jayathilake and Smith [2020] D. I. Jayathilake and T. Smith, “Understanding the role of hydrologic model structures on evapotranspiration-driven sensitivity,” Hydrological Sciences Journal 65, 1474–1489 (2020).
  * Gajamannage, Paffenroth, and Jayasumana [2020] K. Gajamannage, R. Paffenroth, and A. P. Jayasumana, “A patch-based image denoising method using eigenvectors of the geodesics’ gramian matrix,” arXiv preprint arXiv:2010.07769 (2020).
  * Gajamannage and Bollt [2017] K. Gajamannage and E. M. Bollt, “Detecting phase transitions in collective behavior using manifold’s curvature,” Mathematical Biosciences and Engineering 14, 437–453 (2017).
  * Bergström  _et al._ [1995] S. Bergström _et al._ , “The hbv model.” Computer Models of Watershed Hydrology. , 443–476 (1995).
  * Chiew  _et al._ [2002] F. Chiew, M. Peel, A. Western, _et al._ , “Application and testing of the simple rainfall-runoff model simhyd.” Mathematical Models of Small Watershed Hydrology and Applications , 335–367 (2002).