1

Deep Recurrent Neural Networks for TimeSeries Prediction
Sharat C Prasad, Member, IEEE, and Piyush Prasad

Abstract— Ability of deep networks to extract high level
features and of recurrent networks to perform time-series
inference have been studied. In view of universality of one hidden
layer network at approximating functions under weak
constraints, the benefit of multiple layers is to enlarge the space
of dynamical systems approximated or, given the space, reduce
the number of units required for a certain error. Traditionally
shallow networks with manually engineered features are used,
back-propagation extent is limited to one and attempt to choose a
large number of hidden units to satisfy the Markov condition is
made. In case of Markov models, it has been shown that many
systems need to be modeled as higher order. In the present work,
we present deep recurrent networks with longer backpropagation through time extent as a solution to modeling
systems that are high order and to predicting ahead. We study
epileptic seizure suppression electro-stimulator. Extraction of
manually engineered complex features and prediction employing
them has not allowed small low-power implementations as, to
avoid possibility of surgery, extraction of any features that may
be required has to be included. In this solution, a recurrent
neural network performs both feature extraction and prediction.
We prove analytically that adding hidden layers or increasing
backpropagation extent increases the rate of decrease of
approximation error. A Dynamic Programming (DP) training
procedure employing matrix operations is derived. DP and use of
matrix operations makes the procedure efficient particularly
when using data-parallel computing. The simulation studies show
the geometry of the parameter space, that the network learns the
temporal structure, that parameters converge while model output
displays same dynamic behavior as the system and greater than
.99 Average Detection Rate on all real seizure data tried.
Index Terms— Recurrent Neural Networks, Multi-layer
Neural Networks, Artificial Neural Network, Neural Prosthesis

T

I.

INTRODUCTION

IME-SERIES are found in many areas including
communication, health and finance. Large-scale computer
and communication networks generate time-varying metric
that characterize the state of the networks. Audio signal
arriving at a microphone or radio-frequency signal arriving at
a receive antenna are series of values. In the area of health
Electrocardiograph, electroencephalogram, etc all are or
generate time-series data. End-of-day face-values of different
Date first submitted 8/15/2014.
Sharat C Prasad is Best3 Corp, Saratoga, CA 95070 (e-mail:
scprasad@ieee.org).
Piyush Prasad is with the Dept of Biomedical Engineering, University of
Washington, ST Louis, MO 65030 (e-mail:piyushsprasad@gmail.com).

financial instruments are time-series. In all these cases, it is of
interest to discern patterns, detect behavior that differs from
ordinary and predict future behavior.
Many applications have understandably stringent
requirements. We pick implantable epileptic seizure
suppression electro-stimulator driven by automatic seizure
prediction as such an application and to provide the context
for a study of time-series prediction. Epileptic seizure is an
abnormality characterized by brief and episodic neuronal
synchronous discharges with dramatically increased
amplitude.
Electroencephalography (EEG) provides measure of cortical
activity with millisecond temporal resolution. Signal
representations and processing of EEG recordings has been
considered since early days of EEG. Recordings are band-pass
filtered (e.g. with a pass band of 1-70Hz) and digitized (e.g. at
200 samples/second and 12-bit resolution).
The current gold standard for determining onset times is the
trained eye of an encephalographer. Expert encephalographers
look for rhythmic activity. Brain activity during a seizure
displays a great deal of variation from patient to patient and
even from one incident to the next for the same patient.
Implantable neural prosthetic that deliver electrical
stimulation on-demand have emerged as a viable and preferred
treatment for a portion of the population of seizure patients.
Up until recently for reasons of feasibility of implementation,
systems normally take the form of bed-side equipment
communicating with implanted electrodes.
Size is a consideration for any implantable device.
Furthermore both on account of power availability and of
damage to surrounding tissue by the heat produced by the
implanted device, power consumption is a consideration.
Detection algorithms employ characteristics of or features
extracted from EEG signal for detection.
Both feature computation and prediction modules contribute
to the overall size and power requirement. Feature
computation and prediction employ algorithms required to
achieve their objective. Each feature computation may
compute a distinct feature and only this feature.
Savings result if same computation element is reused for
sub-phases of the same computation phase (e.g. different
feature computations). Greater savings result if same
computation element is reused across phases (e.g. feature
computation and prediction).
Given the variability in signature in EEG of seizure even for
the same patient, a rich variety of features has been designed
to assist with the task of seizure detection [17][9]. Given the

2
consideration of surgery, it is desirable that the implant
include support for any feature computation that is or may
become necessary1. Support for comprehensive set of features
is made complicated by the possibility that the optimal set of
features may differ for the same patient from one stage of
treatment to a following stage.
Given the vector of inputs
, vector of features
is
computed. System is modeled as having internal state

and

transition function

and

measurement function

such that
such that

In defining the internal state

, attempt is made to satisfy the
.

As is usual to the type of problem at hand, the parameters of
the prediction algorithm may need to be set by a training
phase. Also, as elsewhere, there are two options to training offline and online. A consideration with offline trained
algorithms is that as parameters change, the offline training
needs to be repeated (or also continuously run in parallel) and
the detection algorithm needs to be configured afresh with the
new parameters.
Systems employing online training can be considered to
continuously incrementally train themselves and therefore able
to adapt.
We will refer to frameworks where given the sequence of

U[

k+1 ]

and

= { u1 , u 2,…, u k+1}

observations

[k]

= { z1 , z 2,…, z k } where each
is
,
is a
known – the measurement - function but
s are not
€
observable,
notation
stands for estimate of
given
Z

€

and

, we first predict

is available, compute

{xˆ

k+1| k+1

and then, when

, xˆ k| k+1,…, xˆ k−q+2| k+1

} as predict-

correct frameworks. The ability of predict-correct frameworks
to “go back” and compute an improved estimate
than

€

computed earlier offers the possibility that

occurrence of an easier to detect event at instant
can be
used to learn correct mapping when
were either not
available or were not informative.
i=M
The problem is one of computing parameters w = { w i } i=1
of a parameterized function

such that when it is applied to

the training input data X = { x j }
from values Zˆ = zˆ j = f w x j

{

( )}

j=N
j=1

j=N
j=1

and supplied output
values Z = { z j }
€

some metric computed
€
returned by the function

j=N
j=1

is optimized.

For example, the problem can be expressed as
€

€
1 Another application
area with this characteristic is Deep Space
Exploration. Once a explorer is launched with certain, preferably learning,
computation engines, these computation engines have to learn to solve the
problem even if the characteristics of the problems evolve as explorer travels
deeper and deeper into outer space.

(1)

Here
is a trade-off constant that trades-off training error
minimization for regularization.
As expressed above the problem is called ridge-regression
and has a simple closed form solution
.

, are defined.

Markov condition

inputs

.

(2)

The above solution represents an off-line training method and,
as introduced above has following assumptions:
•
is linearly dependent on
,
•

s are i.i.d.,

•

components of

•

given training pairs (

are not linearly correlated,
,

) are noise free and

objective and regularization are as in Ridge
Regression.
We know each of these to not hold for Epileptic Seizure
prediction. Adaptability and therefore online training is
desirable in seizure prediction methods. Only a highly
complex feature can hope to achieve linear relation to
prediction output. Individual feature vectors in a sequence of
feature vectors are not independent. Training data is usually
significantly noisy. An objective such as maximizing average
detection rate is not maximizing linear combination of
distances.
When components of
are not completely orthogonal,
methods exist to orthogonalize them [13].
When
s are not i.i.d., 'Time Series' methods consider
that
instant'

is value at 'time instant'

,

is value at 'time

and so on. These methods take covariance

explicitly into account.
Assumption of linear relation of

to

can be

eliminated by employing a Neural Network.
Feedback interconnections within Recurrent Neural
Networks impart them the ability to retain information and
enable them to perform inference tasks such as Epileptic
Seizure Prediction involving time series such as
Electroencephalographic data of brain activity.
Deep Neural Networks can automatically extract relevant
high-level characteristics making it unnecessary to have to
manually engineer characteristic extraction and avoiding loss
of performance from characteristics being not matched to task
at hand. However, most importantly, automatic feature
extraction, addresses the need for adaptability whereby
features suited to a patient are used.
However both, the network being deep and network
incorporating feedback, add to the difficulty of training the
network. Approaches considered for training deep networks
include combining and following unsupervised pre-training
with supervised training [4]. At core of both of these training
algorithms is usually the backpropagation gradient descent.

3
Among the earliest methods considered for training
Recurrent Networks is the Back-Propagation Through Time
(BPTT) method [24] that can be considered the analogue of
the Back-Propagation method for Feed-Forward Networks.
However the BPTT method had been believed to be prone to
getting caught in many shallow local minima on the error
hyper-surface in the parameter hyper-space. Traditionally, for
local minima avoidance, unit activation and/or weight update
are made stochastic and the activation or the weight are
assigned, rather than the exact computed value, a value drawn
from the applicable distribution with the computed values as
moments of the distribution.
Furthermore limitations stem from approximating the nonlinear dynamics of systems as (locally) linear for purposes of
deriving the moments of the posterior state distribution and
from approximating the state distributions as Gaussian.
It is to be noted that greater the non-linearity greater are the
errors introduced by the linear approximation. In case when
process and measurement functions are deep neural networks,
if the transformation from external inputs thru the collection of
layers to final outputs is considered, the transformation
represents a large degree of non-linearity.
Computing the derivatives is tedious and is considered
another shortcoming of methods requiring derivatives.
It is needed to address both the eventual and the speed of
convergence. The model of a system, that when supplied with
some input with stationary statistics, is producing the same
time varying output as the system, cannot be considered
converged if one-step updates are non-zero.
The main contribution of the present work are (i) Deep
Recurrent Neural Network as a unified solution for automatic
extraction of relevant features, automatic construction of
internal state and prediction for time-series data, (ii)
theoretical derivation of the broader space of dynamical
systems enabled by both more than 1 hidden layer and backpropagation through time by greater than 1 instant (iii) the
Back-Propagation Thru Time-and-Space training method as a
solution to hard problem of training deep recurrent networks,
(iv) the dynamic programming formulation of the weight
computation and (v) small, low-power and adaptive epileptic
seizure prediction solution.
In outline, this report presents in Section 2, Deep Recurrent
Neural Networks and their Back-propagation Thru Time and
Space Training. Section 3 discusses universality of Recurrent
Neural Networks and the bound on error in approximation by
RNN. Section 5 presents simulation study of training of and
prediction using a Deep Recurrent Neural Network for
automatically extracted features.

Electroencephalography (EEG) provides measure of cortical
activity with millisecond temporal resolution. It is a record of
electrical potential generated by the cerebral cortex nerve cells
(Figure 1). When epileptic brain activity occurs locally in a
region of the brain and it is seen only in a few channels of
EEG and, when it happens globally in the brain, it is seen in
all channels of EEG.
An example setup for EEG recording [21] may include
Ag/AgCl disc electrodes placed in accordance with one of the
electrode placement systems (e.g. the 10-20 international
electrode placement system [16]). Electrodes may be placed
on the scalp, or may be leads that are implanted in the
cranium. These may be implanted in or on the cortex at the
seizure onset zones. An example is the sensor-neurostimulator device from NeuroPace [23] that includes an
implantable battery-powered neuro-stimulator and quadripolar
depth or strip leads that are implanted in the cranium and are
used for both sensing and stimulation.
Signal representations and processing of EEG recordings
has been considered since early days of EEG. Recordings are
band-pass filtered (e.g. with a pass band of 1-70Hz) and
digitized (e.g. at 200 samples/second and 12-bit resolution).
Representations based on Fourier Transform have been
popular. Four characteristic bands - delta (<4 Hz), theta (4-8
Hz), alpha (8-13 Hz) and beta (13-10 Hz) are identified in the
Fourier Transform.
The current gold standard for determining onset times is the
trained eye of an encephalographer. Expert encephalographers
look for rhythmic activity. Epileptiform discharge may be seen
in the form of spikes repeating at a frequency of 3-Hz and
wave complex in absence seizure (Figure 2). Phase
synchronization and changes in it, are believed to accompany
seizures. Reported measures representative of seizure also
include spectral power, synchronization or amplitude deviates
from normal bounds.

II. TIME-SERIES AND PREDICTION
A. Epileptic Seizure
Epileptic seizure is an abnormality characterized by brief
and episodic neuronal synchronous discharges with
dramatically increased amplitude. Seizures are generated by
abnormal synchronization of neurons which is unforeseeable
for the patients.

Figure 1: Normal brain activity.

Brain activity during a seizure displays a great deal of
variation from patient to patient and even from one incident to
the next for the same patient. A patient may exhibit low
voltage fast gamma activity, may have a beta frequency burst

4
just prior to high frequency activity, may have onset activity in
alpha bands characterized by rhythmic round sinusoidal
waves, may have sharp waves in the delta-theta range
followed by high amplitude spiking activity, may have one of
these onset activities or voltage gamma activity preceded by
spike-and-wave activity or may have activity from the above
entire range preceded by semi rhythmic slow waves [17].
There are two benefits to the Time Series analysis of EEG
data advocated here - (1) offline to determine whether pre-ictal
states exist enabling long-term advance prediction of seizure

While wireless communication between implanted electrode
and external control unit is possible to simplify the implanted
electronics, this may limit fine-grained analysis as the possible
data rates may remain of the order of 100s of Kbps and
therefor too low for some time to come [25]. While new
wireless powering technology is coming into being, the rate of
energy transferred may remain low for some time to come
[18]. Hence size and power consumption need to be
considered along side detection efficiency.
Detection algorithms employ characteristics of or features
extracted from EEG signal for detection. Feature may be
computed on the raw time-domain or frequency or another
transform domain representation. For example, sequence of
samples can be divided into time windows, frequency
transform can be computed for each time window to obtain
spectral power in a few different frequency bands, the vector

Figure 2: Brain activity during seizure.

and (2) online employing information gained from (1) for
preventative therapy in the form of electro-stimulation that
have been seen to be able to suppress or limit the seizure if
applied in a timely manner. Electro-stimulation is believed to
reset brain dynamic from pre-ictal to inter-ictal state.
Implantable neural prosthetic that deliver electrical
stimulation on-demand have emerged as a viable and preferred
treatment for a portion of the population of seizure patients.
While timely application of electro-stimulation is seen to
suppress seizures, extended or more than necessary electrostimulation is seen to damage neuronal tissue and have such
adverse consequences as requiring increasing amounts of
stimulation to have effect and stimulation induced suppression
of neuronal excitability. Up until recently for reasons of
feasibility of implementation, systems normally take the form
of bed-side equipment communicating with implanted
electrodes. In comparison self-contained neural prosthetic will
not require the patient to be bed-bound. Responsive neurostimulation, by virtue of being timely and limiting the amount,
may also significantly increase the effectiveness of therapy.
The advent of very low-power battery-operated electronics has
made possible implantable electronics with closed-loop
working of detection and stimulation.
Traditionally detection algorithms have emphasized
detection efficacy and have disregarded power consumption
and size. Size is a consideration for any implantable device.
Furthermore both on account of power availability and of
damage to surrounding tissue by the heat produced by the
implanted device, power consumption is a consideration.

Figure 3: Organization of detection sub-system.

of spectral powers can be a feature vector and a sequence of
such vectors may be analyzed by a detection algorithm. Figure
3 shows the overall organization of the detection sub-system
comprising of feature computation and prediction modules.
Both feature computation and prediction modules contribute
to the overall size and power requirement. Feature
computation and prediction employ algorithms required to
achieve their objective. Each feature computation may
compute a distinct feature and only this feature.
In above, we gave examples of characteristics expert
encephalographers look for in the EEG. Some of these
characteristics were described in terms of frequency domain
representation. To make use of such characteristic, frequency
domain representation needs to be computed. Efficient
frequency transform computation requires specialized
algorithms and space, time and power. If realized as circuits,
the circuits may not serve to compute a different
representation.
Time-domain features avoid the cost of transform
computation but require specialized algorithms.
One noteworthy aspect of the transform computation

5
algorithms, that we will use to argue in favor of our thesis, is
that their regular structure enables easy to visualize trading-off
of space (number of computation engines, amount of memory)
for time and remain acceptable as long overall latency remains
within limits. To save space and power, computation may be
performed iteratively, reusing the same computation element.
Savings result if same computation element is reused for subphases of the same computation phase (e.g. different feature
computations). Greater savings result if same computation
element is reused across phases (e.g. feature computation and
prediction).
Given the variability in signature in EEG of seizure even for
the same patient, a rich variety of features has been designed
to assist with the task of seizure detection [17][9]. There have
been studies of efficacy of different detection features. Given
the consideration of surgery, it is desirable that the implant
include support for any feature computation that is or may
become necessary2. Support for comprehensive set of features
is made complicated by the possibility that the optimal set of
features may differ for the same patient from one stage of
treatment to a following stage.
Advent of small low-power computer-on-chip [20] makes
them one option to address the feature computation
requirement. However, as numerous studies have shown, size
and power requirement of general purpose programmable
solutions always exceed custom optimized solutions.
There have also been studies of attempts [9] to improve
performance by using detection features in combination.
Latter has shown that when all of a set of detection features
are used best performance results. However the benefit
€ from
using more than a certain number of features may not justify
their incremental cost when size and power constraints exist
[17].
B. Time-series prediction
Given the vector of inputs

, vector of features

is

computed. System is modeled as having internal state

and

transition function

and

measurement function

such that
such that

In defining the internal state

, are defined.

, attempt is made to satisfy the

Markov condition

.

As is usual to the type of problem at hand, the parameters of
the prediction algorithm may need to be set by a training
phase. Also, as elsewhere, there are two options to training offline and online. An offline (batch) training algorithm is
given feature vector

and target

pairs

{( x , z )}
i

i

i= N

.

i=1

The training goes over all the data to determine detection
parameters. A consideration with offline trained algorithms is
that as parameters change, the offline training needs to be

€
2 Another application area with this characteristic is Deep Space
Exploration. Once a explorer is launched with certain, preferably learning,
computation engines, these computation engines have to learn to solve the
problem even if the characteristics of the problems evolve as explorer travels
deeper and deeper into outer space.

repeated (or also continuously run in parallel) and the
detection algorithm needs to be configured afresh with the
new parameters.
On the other hand, an online (incremental) algorithm
initializes its parameters to certain initial values, is given
either just one feature vector
or one feature-vector-target
pair

. An online algorithm only produces a prediction

if only given a feature vector or also updates its parameters
if given a feature-vector-target pair

. To start with, an

on-line trained algorithm may produce poor predictions. As
the data characteristics are learnt, predictions improve.
Systems employing online training can be considered to
continuously incrementally train themselves and therefore able
to adapt. One issue with online training is that for supervised
training, target
corresponding to each feature-vector
is
required. A second issue is that, on one hand, recent featurevector-target pairs may cause invariant characteristics to be
over-written by transient characteristics and on the other hand,
recent pairs may not cause transient characteristics to be
overwritten by new transient characteristics.
We will refer to frameworks where given the sequence of
inputs

U [ k+1] = {u1,u2 ,…,uk+1}

and

Z [ k ] = {z1, z 2 ,…, z k } where each

observations

is

,

is a known – the measurement - function but

€observable, notation
and

s are not

stands for estimate of

, we first predict

is available, compute

{xˆ

given

and then, when

}

, xˆ k| k+1,…, xˆ k −q+2| k+1

k+1| k+1

as

predict-correct frameworks. The ability of predict-correct
frameworks to “go back” and compute an improved estimate
than
computed earlier offers the

€ occurrence of an easier to detect event at
possibility that
instant
can be used to learn correct mapping when
were either not available or were not informative.
Metrics called specificity (
Average Detection Rate (

(

), sensitivity (
) and
) are used. Sensitivity is

)

defined as SEN ≡ Y + / Y + + N − where
true positives and

is the count of

is the count of false negatives.

Specificity is defined as
€
the count of true negatives and
positives. Finally

where

is

is the count of false
.

The problem is one of computing parameters
of a parameterized function

such that when it is applied to

the training input data U = { u j }
from values Zˆ = zˆ j = f w x j

{

( )}

j=N

€

returned by the function

j=1

and supplied output
values Z = { z j }
€

€

some metric computed

j=1

j=N

j=N
j=1

is optimized.

€

€

€

€

6
For example, it may be desired that sum of squared errors is
minimized. Then the problem is, find
where
(3)
Having computed

prediction output. Individual feature vectors in a sequence of
feature vectors are not independent. Training data is usually
significantly noisy. Objectives such as maximizing average
detection rate are not maximizing linear combination of
distances.

, prediction is made as

.
For the problem under consideration,
, and a
problem encountered is that of over-fitting wherein a
is
determined that well minimizes the sum of squared errors for
the training data (by learning too well the characteristics
present) but then performs poorly on other valid data (that
may exhibit a broader set of characteristics). To counter this, a
regularization term may be added to the minimization and the
problem may be written as

argmin ⎛
ˆ=
w
∑⎜⎝ f w u j − z j
w

( ( )

) + λg(w)⎞⎟⎠
2

If it is assumed that the relation of output
input is linear, i.e.

(4)

to corresponding

m

z j = ∑ uij w i = uTj w j ,

(5)

i=1

s are drawn from independent identical distributions (
does not depend on any

,

and

for all

drawn from the same distribution), components of

, are
are not

Figure 4: Prediction with manually engineered features

linearly correlated and lastly the regularization function
simply takes the sum of square of individual elements, i.e.
, the problem can be expressed as
2
argmin ⎛
⎞
ˆ=
w
⎜ f w u j − z j + λg( w)⎟
∑
⎝
⎠
w
argmin
T
=
(Uw − z) (Uw − z) + λwT w .
w

( ( ) )

(

)

(6)

Here
is a trade-off constant that trades-off training error
minimization for regularization.
As expressed above the problem is called ridge-regression
and has a simple closed form solution
.

(7)

The above solution represents an off-line training method and,
as introduced above has following assumptions:
•
is linearly dependent on
,
•

s are i.i.d.,

•

components of

•

given training pairs (

are not linearly correlated,
,

) are noise free and

objective and regularization are as in Ridge
Regression.
We know each of these to not hold for Epileptic Seizure
prediction. Adaptability and therefore online training is
desirable in seizure prediction methods. Only a highly
complex feature can hope to achieve linear relation to
•

Figure 5: Metrics SPC, SEN and ADR versus illustrative manually
engineered features on a sample of real data shows SPC, SEN and ADR with
respect to four manually engineered features [17] – Energy, Hjorth Variance,
Coastline and RMS Amplitude.

Complex objectives rarely have solution in closed form. As
long as the objective is differentiable, we can employ a
member of the Gradient Descent family of algorithms for
solution.
Convex
functions,
characterized
by
are able to express a wider set
of objectives. Carefully designed members of the Gradient
Descent family can optimize convex objectives, can
accommodate samples that are not independent [7] and even
overcome noise and certain forms of non-convexity by taking
advantage of the structure in the objective [15].

7
When components of

are not completely orthogonal,

methods exist to orthogonalize them [13]. Most regression
packages detect when some components are linearly correlated
and either suggest components to eliminate or set respective
coefficients to zero.
Ignoring for the moment that
s are not i.i.d. and
employing manually engineered features, we can use
organization shown in Figure 3 to perform prediction. The
organization uses a single neuron. The weights of the neuron
can be trained using a member of the Gradient Descent family
of algorithms. This organization does not compute any further
higher layer features beyond those manually engineered and
makes predictions that are non-linear but not highly non-linear
functions of the features.
Figure 5 shows SPC, SEN and ADR with respect to four
manually engineered features 0 – Energy, Hjorth Variance,
Coastline and RMS Amplitude. The non-smooth and gradientdiscontinuous nature of the plot is from, among others,
presence of noise and use of discrete 'counts' of 'events' – a 0 if
did not occur and a 1 if it did - in the definition of the SPC and
SEN.
When
s are not i.i.d., 'Time Series' methods consider
that
instant'

is value at 'time instant'

,

is value at 'time

and so on. These methods take covariance

explicitly into account.
Assumption of linear relation of

to

can be

eliminated by employing a Neural Network. Often inference
requires computing an unknown and, possibly highly, nonlinear function of some characteristics of the data. Neural
Networks can efficiently find nonlinear mappings.
Feedback interconnections within Recurrent Neural
Networks impart them the ability to retain information and
enable them to perform inference tasks such as Epileptic
Seizure Prediction involving time series such as
Electroencephalographic data of brain activity.
Deep Neural Networks can automatically extract relevant
high-level characteristics making it unnecessary to have to
manually engineer characteristic extraction and avoiding loss
of performance from characteristics being not matched to task
at hand. However, most importantly, automatic feature
extraction, addresses the need for adaptability whereby
features suited to a patient are used. Furthermore, as the
treatment progresses, the features used can also evolve to suit.
However both, the network being deep and network
incorporating feedback, add to the difficulty of training the
network. Approaches considered for training deep networks
include combining and following unsupervised pre-training
with supervised training [4]. The unsupervised pre-training
tries to minimize energy as low-energy states correspond to
parameters attaining values that capture structures present in
data. The following training then is intended to further evolve
the parameters to capture the desired input-output mapping. At
core of both of these training algorithms is usually the
backpropagation gradient descent.

Figure 6: Automatic High-Level Feature Extraction Using Deep Recurrent
Neural Network

Among the earliest methods considered for training
Recurrent Networks is the Back-Propagation Through Time
method [24] that can be considered the analogue of the Back
Propagation method for Feed-Forward Networks. However the
BPTT method had been believed to be prone to getting caught
in many shallow local minima on the error hyper-surface in
the parameter hyper-space. Traditionally, for local minima
avoidance, unit activation and/or weight update are made
stochastic and the activation or the weight are assigned, rather
than the exact computed value, a value drawn from the
applicable distribution with the computed values as the
moments of the distribution.
It was first noted in late 80s that the Neural Networks are
not unlike non-linear physical systems and that there exists a
well-developed theory of estimating hidden state and
estimating parameters of linear systems from noisy
observations and somewhat less developed similar theories for
nonlinear systems. A number of methods adhering to PredictCorrect framework (e.g. the celebrated Kalman Filter) have
been developed over time that are optimal in the sense of
Minimum Mean Squared Error for state only, parameter only
and joint state and parameter estimation problems for linear
systems. Optimal estimation is no longer computationally
tractable when the system is nonlinear.
The realization that Recurrent Neural Networks are nonlinear systems led to use of Extended Kalman Filter, popularly
used for non-linear system identification, parameter
estimation, etc for training Recurrent Networks. The Extended
Kalman Filter has a set of well known limitations. The
limitations stem from approximating the non-linear dynamics
of systems as (locally) linear for purposes of deriving the
moments of the posterior state distribution and from
approximating the state distributions as Gaussian.

€

8
Time-series prediction has been addressed by Statistics,
Dynamical Systems, Machine Learning, Signal Processing
and Automatic Control communities. The need to address
larger and more complex systems has led beyond BackPropagation Through Time and first order term (in Taylor
expansion) only Extended Kalman Filter [11] to including
second order terms [14], avoiding explicit use of Jacobian and
Hessian [14], Quadrature [1] and Cubature [1] forms, the
Unscented [12] form, the marginalized Sigma Points [19]
form, the Particle Bayesian [22] form, the Variational
Bayesian [22] form and the Gaussian Process [3] form. While
an improvement results, a conceptual and sometimes also a
computational complexity cost has to be paid.
It is to be noted that greater the non-linearity greater are the
errors introduced by the linear approximation. In case when
process and measurement functions are deep neural networks,
if the transformation from external inputs thru the collection of
layers to final outputs is considered, the transformation
represents a large degree of non-linearity.
Computing the derivatives is tedious and is considered
another shortcoming of methods requiring derivatives.
It is needed to address both the eventual and the speed of
convergence. The model of a system, that when supplied with
some input with stationary statistics, is producing the same
time varying output as the system, cannot be considered
converged if one-step updates are non-zero.
C. One-step and long-term time-series prediction
The problem comprises of being given an input time-series

instead
be
approximations

obtained

as

a

sequence

of

,
,
, etc.
In the context of training multi-layer recurrent neural
networks, where weight update is computed as
, this fact can be taken advantage of by
including in the state, in addition to prior outputs of hidden
layers, the current output of visible and hidden layers.
E. High order Markov Processes
Given a process that receives input
produces output

at instant

and

, the process may employ internal state

to remember as much of past history as it needs to so that
it can compute output

considering only

and

.

Expressed another way,
and system output is considered independent of past input
given the current internal state. Such processes are called firstorder. However this may require, not unlike carefully
manually engineering features, carefully manually engineering
the internal state that system must remember. A solution is to
require system to remember history going back
instants.
Then
and
the process is called of order

.

and noisy observations of corresponding output
time-series

where observations may not be

available at each instant of time
and the required
result being
. As an example, the input time-series
may be intra-cranial electrode sensed electric potential, the
output time-series values may only be available for certain ,
may be whether the patient was experiencing an epileptic
seizure3 and the required result is whether the patient is about
to experience an epileptic seizure instants from now.
D. Extrapolating highly non-linear functions
Consider a function
a

sequence

of

y = f ( x ) = e a1 e

functions:

a2x+b2

+b1 specified as

y = f1 ( y1 ) = e a1 y1 +b1

and

a 2 x+b 2

.
y1 = f 2 ( x ) = e
We are given
€ value

of

at

and it is desired to

extrapolate value at €
.
A linear approximation can be done and we can write
.
It is easy to show that a more accurate approximation can
3 Note that in an offline training situation, the EEG data may be
examined by an expert encephalographer and intervals corresponding to
seizure determined. In an online training context, the parameters may be
initialized to values determined offline and predictions used only after
stability is reached.

III. DEEP RECURRENT NEURAL NETWORKS
Neural Networks comprise of one or more inputs, one or
more neurons (referred to in this work as units), one or more
outputs, directed or undirected interconnections among units, a
weight associated with each undirected connection of or each
directed connection into a unit (from an input or another unit)
and each unit computing its activation by applying its
activation function to weighted sum of inputs.
Feed-Forward Neural Networks are arranged in layers with
units in a layer receiving inputs from only units in previous
layer (or outside in case of the input layer). Feed-Forward
networks process data presented at an instant independent of
data presented at previous instants and so are unable to
perform inference on time series data. In a Recurrent Network,
some of the interconnections are from units in a layer through
a delay to units in the same or an earlier layer.
As both stacking of recurrent layers and back-propagating
have the effect of extracting information from input or from
information extracted from input to consider alongside input,
it is a question as to whether (i) either is more effective than
the other and (ii) having both together adds to the
approximation ability of the network.
Each unit in a layer associates a weight to and has as inputs
a fixed bias, output of all units in the previous layer (external
inputs in case of the visible layer) and, in case of the recurrent
layer, its own outputs delayed by a time instant. Output of a

€

€

€

€

9
unit

is the result of applying a nonlinear function

weighted sum of inputs

:

to the

.

The network has a visible layer,

hidden layers and an

output layer of units (neurons) with, respectively,
and
inputs

units in the layer. The network receives

,

and produces outputs

Furthermore,

is the

the visible layer,
the unit

is the input

,

) in
is

is the weight from unit

€,

is

(xˆ ( ) ) . Final outputs of
k

i

the network are also considered collected into a vector. This
vector at instant is denoted as
.
A significant component of any method making use of the
€ computations
derivatives are the derivative

( )
∂( xˆ ( ) )
∂( zˆ ( ) )
( )
(H ) = ∂ xˆ ( )
( )
(F )
( k)

ij

=

∂ xˆ ( k)

i

k

(12)

j

k

ij

i

k

(13)

j

to unit

of the node

l, j

of the same is denoted as

k

.

is the vector with elements

y h ( − k) is the output at (past) instant

hidden layer

.

in the hidden layer ,

in the output layer,

the output of unit

,

-th unit (or alternatively unit

is the unit

As a convention,
,

,

element

€

. Notionally, associated with the unit

in
in

hidden layer
are
memories that retain the previous
outputs y ( − B ) through y ( −1) of the unit.
is referred to as

€

h l, j

€

h l, j

the extent of backpropagation.
Inputs to each unit in any layer include a fixed bias 1. We
use the same symbol as used to denote the output signal vector
of a layer, €
but with a tilde accent, to denote a signal vector
with 1 (corresponding to the fixed bias) as an additional and
the first element and remaining elements same as the original
signal vector. Then
y v = Wv,u y˜ u
(8)

yh

NL

= WhNL , v y˜ v + Wh

NL

( −1)

, h NL

( −1)

h NL

NL

(9)

NL

y h ( −1) = Wh ( −1) ,v ( −1) y˜ u( −1) + W
NL

y h ( −1)
,hN( −2)
L

y h ( −2)

(10)

NL

Reader may notice that we have chosen to qualify, rather
than the symbol denoting the output signal of a node, the
symbol denoting the node with the superscripts
denoting
value of “the node at past instant ” rather than “value at past
instant ” of the node. This is deliberate and is to suggest that
we consider the network to have been unfolded back and there
to exist
copies, one corresponding to each of
instants
spanning the extent of the backpropagation, of the complete
“sub-network” comprising of the visible and each of the
hidden layers. As a result, the notation allows not only vectors
to be distinguished from the vector
, but
also weight matrices
distinguished from

Wv ( −1) , u( −1) ,, Wv ( − B ) , u( − B )

to be

. If this distinction was not made (10)

above will be

y h ( −1) = Wh ( −1) , v ( −1) y˜ v ( −1) + Wh ( −1) , h ( −2) y h ( −2)
NL
N€
NL
NL
NL
L

(11)

denoted as

and

Outputs of all the units as well as the weights associated
with the inputs to the units are considered collected into a
joint-state-and-parameter-vector. This vector at instant
is
, estimate of the vector is denoted as

Figure 7: Back-Propagation Through Time and Space in Deep Recurrent
Neural Networks

Visualization of dependencies considered in the derivative
computation is eased by considering the nodes making up the
network to be arranged as columns of stacks of layers of linear
array of nodes (Figure 7). The input layer is at the bottom, the
hidden layers in the middle and the output layer at the top. The
rightmost column corresponds to current instant and the
columns to left to successively previous instants. A notion of
distance can be associated with any pair of units and of a unit
and a weight. A unit in the (present-instant, output-layer) is at
a distance of 1 from any unit in the (present-instant, hidden
layer), at a distance 2 from any unit in (present instant, visible
layer), also at a distance 2 from (previous instant, hidden
layer) and so on. Similarly the instant
weight
at a distance 1 from (instant
, hidden layer

) unit

is
,

10
distance 2 from any unit in the (present instant, hidden layer
), also distance 2 from any unit in the (instant
, layer
), distance 3 from any unit in the (present instant, layer
), also distance 3 from any unit in (instant

, layer

), etc.
Output at instant
depends on parameters (weights)
€ and
state at . The state at instant
in turn depends on weights
and state at
. Once weights have been updated at instant €
, a new output at
can be computed. To compute this
output, the state at
is required. However the state at
itself can be recomputed using the updated weights and
€
the state at
and so on.
Given the error in output at , part of the error can be
attributed to errors in weights and part to error in
€ the
(estimated) state at . The error in state at
in turn depends
€
in part on errors in weights and part on error in state at
.
If the mapping from input to output is invariant, the weights
at instants
and
are the same. However if the
mapping itself is evolving then the two sets of weights are not
the same.
€
Once the weights corresponding to present and to each past
instant are made different and the present error is computed,
weights at instant
can be overwritten with weights
€ at
instant
and present instant weights can be updated to next
€
instant weights. Alternatively, the back-propagation itself can
be used to compute weights at present as well as each retained
past instant.
The mapping

may be time-invariant,

may need to be learnt and, once learnt, does not change.
€
Alternatively the mapping may continuously evolve.
€
Computing only next instant weights by updating present
instant weights and overwriting instant
weights with
instant
weights is consistent with models that allow
€ the
mapping to be learnt and then assume them to remain
constant. Computing all the weights using back-propagation is
€
consistent with models where mapping may continuously
evolve.
A. Universality of single hidden layer feed-forward network,
bounds on approximation error and multilayer networks

depends on the norm of the Fourier Transform of the function
being approximated
2
2
c f ' = 2rC f = 2r ∫ |ω|| F˜ ( dω )|
(14)

(

(

)

)

where F˜ ( w) Fourier distribution of function
,
i
θ
w
)
such that f ( x) = e iwx F˜ ( dw) , F˜ ( dw) = e ( F ( dw) ,

,
is

∫

the magnitude distribution,

the phase at frequency

and w the vector of frequencies (ω1 ,, ω d ) .
€ Let
€ of functions
be the class
on
for which
C f = ∫ |w|| F˜ ( dw)| < ∞ and
the class of functions
on
€
d
R for which C f = ∫ |w|| F˜ ( dw)| < C .
Then, for every choice of fixed basis functions
[2]

()

€
C ⎛ 1 ⎞
sup s( f,span( h1, h2 ,, hn )) ≥ κ ⎜ ⎟
d ⎝ n ⎠
f ∈ ΓC
where
s( f,G ) = inf s( f, g)

1
d

(15)
(16)

g ∈G

s( f, g) =

2

∫ ( f ( x) − g( x)) dx

(17)

Let,
be a bounded set in R d which contains
the set of functions
on R d
f ( x ) = f ( 0) + e iw.x − 1 F˜ ( dw) holds for

∫(

€
iθ w
F˜ ( dw) = e ( ) F ( dw) ,

)

,
for which
and some

€

the set of functions

Rd

on

for which
(18)

∫ |ω| F ( dω) ≤ C > 0
B

where
|w| B = sup{ w.x} ,

€

(19)

x∈B

then for every function
measure

, for every probability

, for every sigmoidal function

, for every

there exists a linear combination of functions

,

of the form

n

f n ( x) = ∑ c kφ ( a k x +bk ) +c 0
k=1

such that [2]

2
2
Single hidden layer feed-forward networks using the
f n ( x ) − f ( x)) µ ( dx) < c f ' / n = 2rC f / n
(20)
(
∫
monotone cosine squashers are capable of embedding as a
€
special case a Fourier network that yields a Fourier series where Br = x ∈ R d | x| ≤ r .
approximation to a given function at its outputs [8] and thus
possesses all the approximation properties of Fourier €series
Hence adding a sigmoidal layer improves the rate from
representations.
1/ d
to
. Will adding another layer improve further?
A feed-forward network with a single hidden layer of € (1/ n)
continuous sigmoidal nonlinear squashers (in addition to linear By how much?
In a one linear layer network we take a weighted sum of
visible and output layers) are complete in the space of
inputs
and add in a bias
where for
functions on
that have continuous and uniformly formed
€
derivatives up through order
[10].
ease of explanation we number the elements of
and
For a network with a hidden layer with
units, the starting at .
approximation error is bounded by
[2] where
Equation (7) bounds below the approximation error of one

(

{

}

)

11
linear layer network.
In a 1 hidden layer network, we take
distinct copies of
one linear layer networks, feed the output of each into a
nonlinearity, take a weighted sum of the outputs of the
nonlinearities and add in a bias

) (
((
⎛
⎞
= b .⎜( f ( x )),, ( f ( x))⎟ +b .
⎝
⎠

f n1 ( x) = bT . φ a1T x + a10 ,, φ a Tn x + a n0

))

T

n 0 ,1

n 0 ,n

(21)

0

Equation (8) bounds below the rate of error. Note that here
we are concerned with inherent representability using multi€ in
layer neural networks, are not concerned with error
determination of weights and assume that values of weights
minimizing error, if they exist, can be found.
€

€

€

2

∫ ( f (x) − f (x)) µ (dx) .

Consider

f n 2 ( x) is the function

n2

computed by the two hidden layers network. The two hidden
€
layers network has (logically)
instances of the one hidden
€
€
€ network as sub-networks.
layer
Let
be the function
computed by the two hidden layers network when it achieves
minimum error and f˜ ( x),, f˜ ( x) be the functions

(

n1 ,1

n1 , n

)

needed to be computed by its one hidden layer sub-networks
when the two hidden layers network achieves minimum error.

(

€

€

€

€

€
€

)

Let f n1 ,1 ( x),,
€ f n1 , n ( x) be the actual functions computed by
€
the one hidden layer sub-networks. Similarly let
f˜n 0 ,1 ( x),, f˜n 0 , n ( x) be the functions needed to be computed
€
by their one linear layer sub-networks when the one hidden
layer network achieves minimum error relative to
f˜ ( x),, f˜ ( x) and f ( x),, f ( x) be the actual

(

(

n1 ,1

n1 , n

)

(

n 0 ,1

)

n0 , n

functions computed by the one linear layer sub-networks.
Let us denote
⎛
ψ n i ( x) = ⎜φ f n i ,1€
( x) ,, φ f n i , n ( x) ⎞⎟⎠
⎝
€
⎛
⎞
ψ˜ n i ( x) = ⎜φ f˜n i ,1 ( x) ,, φ f˜n i , n ( x) ⎟
⎝
⎠
€
⎛⎜
⎞⎟
+
+
+
ψ n i ( x) = φ f n i ,1 ( x) ,, φ f n i , n ( x)
⎝
⎠
€
⎛⎜
⎞
−
−
−
ψ n i ( x) = φ f n i ,1 ( x) ,, φ f n i , n ( x) ⎟
⎝
⎠
€
Then
€ (22)
f ∗2 ( x) = ψ˜ 1 ( x).c +c 0 ,

(
(
(
(

)
)
)
)

(
(
(
(

)
)
)
)

n

€

f n1 ( x) = ψ n 0 ( x).b+b0 ,

€

∫ ( f ( x) − f˜ ( x)) µ ( dx) < c ' / n

€ (23)

2

€

)

2

/n.

What if the errors of functions computed by the
2
€
2rC f / n
subnetworks
are
smaller
e.g.
versus
€ (1/ d )
κ (C / d )(1/ n) ?

)

Note function computed is weighted sum of
€
f 1 ( x) = ψ n 0 ( x).b+b0

terms
(26)

n

For squared minimum error
⎛⎜ ∗
⎞2
f 1 ( x ) − f ( x ) ⎟ =
inf
ψ n 0 ( x).b+b0 − f ( x)
⎝ n
⎠
f 0 ,, f 0 ,b ,b0
n ,1
n ,n

(

(

)

= ψ˜ n 0 ( x ).b+b0 − f ( x)

)

2

2

(27)

While there can be cases of small2 input errors and errors
adding, such that a large output error results. However in
expectation and averaged over large number of points, smaller
input (sub-network) errors will result in small output
(network) errors.
If
2
∫ ψ˜ 0 (x).b+b0 − f ( x) dµ ( x)

(

=

)

n

∫ (ψ (x).b' +b'

inf
,, fn0,n ,b' ,b0'

n0

f 0
n ,1

0

)

2

(28)

− f ( x ) dµ ( x )

and

)

n

€

(

bounded above by 2rC f

n2

(

+b0

T

T

2

∫ ( f (x) − f ( x)) φ (dx)

subnetworks manage to achieve error

1

n ,k

1

f

n ,k

and

d

(

= 2rC f

)

2

/nd

€ (24)
€

()
1

C ⎛ 1 ⎞ d
(25)
sup s f˜n 0 , f n 0 ≥ κ ⎜ ⎟ .
d ⎝ n ⎠
f ∈ΓC
Hence in spite of the one linear layer subnetworks having

(

)

error bounded below by

, 1 hidden layer

(29)
then
2

∫ (ψ (x) .b+b − f (x)) dµ ( x)
< ∫ (ψ ( x ) .b+b − f ( x )) dµ ( x )
−
n0

0

€
€

(30)

0

Consider
2

2

(ψ (x) .b+b − f (x)) − (ψ ( x) .b+b − f ( x))
= (ψ ( x ) .b+b − f ( x )) − (ψ˜ ( x) .b+b − f ( x))
−(ψ ( x) .b+b − f ( x)) + (ψ˜ ( x) .b+b − f ( x))
= ((ψ ( x) +ψ˜ ( x )) .b+ 2b − 2 f ( x)) ×((ψ ( x) − ψ˜
−((ψ ( x) +ψ˜ ( x)) .b+ 2b − 2 f ( x )) ×((ψ ( x) − ψ˜
+
n0

−
n0

0

0

2

+
n0

2

n0

0

0

2

−
n0

2

n0

0

+
n0

n0

−
n0

n0

(31)

0

0

+
n0

n0

( x)) .b)

0

−
n0

n0

( x)) .b)

>0
(32)
€ each (inner) term of both factors
The last step follows from
of the first (outer) term of€the above expression being greater
than or equal to the corresponding (inner) term of the second
(outer) term.
Hence (11) follows.
Intuitively
difference
⎛⎜
⎞
∫ span⎝φ f n L −1 ,1 (x) ,, φ f n L −1 , n (x) ⎟⎠ dµ (x)

(

∫ span( f

€

2

+
n0

n L −1 ,1

) (

( x),, f n

L −1

,n

between

)

( x)) dµ ( x)

over

and
the

entire

12
hypersphere
function

( 0,1)
f ( x) of

d

can be understood by considering a
which remains equal to some value with

large magnitude
nearly everywhere but dips steeply down
to € at some number of (related to ) points. A linear
€function can only descend and then ascend back up at rate
determined by the linear coefficients and will produce large
error. Note nonlinearity such as
is not convex.
However
convex/concave

is. Furthermore, nonlinearity may still be
between

−( n − 2) π < x < −( n − 1) π

−( n −1) π < x < −nπ .

€

€

and

(vector of weights for the hidden layer unit at instant
also
, .. ,
and we compute
, outputs of hidden layer unit
through
but also

at instants

. When the weights are updated, not only
, .. ,

y H,i,k−1 are recomputed.
instant
history.

) but
, .. ,

are updated and
and

, .. ,

, to be computed at

now make use of recomputed weights and

€
IV. BACK-PROPAGATION THRU TIME AND SPACE
For Back-Propagation Thru Time and Space Gradient
Descent Training, weight
is updated according to
where

is the error,

is the

output of the single output network at the instant when input
has been applied at its input and
is a training inputoutput pair.
Visualization of dependencies considered in the derivative
computation is eased by considering the nodes making up the
neural network to be arranged as columns of stacks of layers
of linear array of nodes (Figure 7). The output layer is at the
top, the hidden layers in the middle and the input layer at the
bottom. The rightmost column corresponds to current time and
the columns to left to successively previous times.
A notion of distance can be associated with any pair of units
and any pair of a unit and a weight. A unit in the (presentinstant, output-layer) is at a distance of 1 from any unit in the
(present-instant, hidden layer
), at a distance 2 from any
unit in (present instant, visible layer), also at a distance 2 from
(previous instant, hidden layer
) and so on. Similarly the
present instant weight

is at a distance 1 from (present

instant, visible layer) unit

, distance 2 from any unit in the

(present instant, hidden layer

Figure 8: Minimum error (over all output layer weights) versus hidden layer
weights for input weights (-1.5, -1.5)

The ability of predict-correct frameworks to “go back” and
compute an improved estimate
than the
offers the possibility that occurrence of an easier to detect
event at instant
can be used to learn a mapping that
produces a better
than when
alone was
available. This is especially relevant to seizure onset
prediction as determining that a seizure is about to occur is
harder for both humans and computers than determining that a
seizure is in progress.

), distance 3 from any unit in

the (present instant, hidden layer

), distance 3 from any

unit in (previous instant, hidden layer

), etc. Hence the

distances here are partially over space and partially over time.
Details of computation of the derivatives are in the
Appendix. These suggest a dynamic programming method for
computing the derivatives comprising of traversing the graph
from inputs to outputs, computing derivative of unit outputs
wrt distance 1 away unit outputs, distance 1 away weights and
distance 2 away weights. Derivatives wrt unit outputs at
distance greater than 1 away and weights at distance greater
than 2 away are computed by taking weighted sums over
derivatives of distance 1 away unit outputs.
A. Adaptation and Autonomous Self Improvement
BPTT implies that at instant

, we not only have

Figure 9: Minimum error (over all output layer weights) versus hidden layer
weights for input weights (-1.5, 1.5)

13
V. EXPERIMENTAL RESULTS
The first set of experiments explored the geometry of the
parameter space first with (hidden layers 1, back-propagation
extent 1), (hidden layers 1, back-propagation extent 2) and
finally with (hidden layers 2, back-propagation extent 1). The
visible, hidden and output layers each had one unit. Units in
visible and output layers had one input from outside or the
previous layer and one bias input. The unit in the hidden layer
had one input from the visible layer, one bias input and one
input that was its own output delayed by one instant. There
was one input to and one output from the network. Hence
there are in total 7 weights in the first, 12 in the second and 10
in the last case. Each weight corresponded to one dimension.
A point in the first parameter space was specified by 7
coordinates. Each coordinate was restricted to one of a set of
discrete values. Range spanned by all of the coordinates was
restricted to a subrange of [-2, 2). For different coordinates,
different number of discrete values evenly distributed across
the range were visited during the systematic exploration.

Figure 12: Error
versus
layer
weights
hidden
and input layer
Figure
11: output
Predicted
versus
noisewith
corrupted
actual
weights corresponding to minimum error

Figure 13: Predicted versus noise corrupted actual
Figure 10: Minimum error (over all output layer weights) versus hidden layer
weights for input weights (1.5, -1.5)

For the first two coordinates, corresponding to the weights
associated with the visible layer unit, typically a small number
(e.g. 2), for the next 3 coordinates corresponding to the hidden
layer, a medium number (e.g. 8) and finally for the last two
coordinates corresponding to the output layer a larger (e.g 32)
discrete values were visited. Example number of values
mentioned above resulted in
different points
in the parameter space to be visited.
For each of the 64 vectors of values of two of the three
hidden layers weights – corresponding to the recurrent input
and input from visible layer - minimum error was determined
separately for each of 4 values of the weight pair
corresponding to weights associated with the visible layer unit.
These are plotted in Figures 8 through 11.
Next with the vector of visible and the hidden layer weights
set to the vector corresponding to minimum error determined
above, the error at the 1024 output layer weight pairs is plotted
in Figure 13.

Figures 8 through 11 have a bearing on the question as to
whether with visible and hidden layer weights initialized
randomly and output layer weights learnt using gradientdescent-backpropagation-through-time suffices to minimize
error to arbitrarily small values. These present empirical
evidence that it is not so. Arbitrary initialization of 'lower'
layers and only learning weights of higher layers lower bounds
the error, such that smaller error cannot be achieved.
Figure 12 has a bearing on the question whether with 'lower'
layer weights set correctly, if the values of the error
corresponding to different values of output layer weights are
considered to lie on a manifold, what does this manifold look
like. As the output layer comprises of sigmoidal units as well,
the manifold is non-convex.
Simulation study was conducted with synthetic and real
data. The synthetic data was the Sine and Mackey Glass
function generated series. In case of the sine function
generated series, the run was conducted upto
.
However during
, prediction was performed but
measurement update was not performed. We see that (Figure
11) predictions match the true data well. Even in case of the

14
Mackey Glass function generated series, which is not simply
periodic, predictions match the true data (Figure 14).
At each instant, input data is received, output is estimated
and, if actual output is made available, the parameters are
updated. It can be argued that with input with stationary
statistics, if parameters do not attain stable values, correct
output being estimated does not mean that input to output
mapping has been learnt. As the plot of evolution of weights
(Figure 14) shows, stable values are attained and the inputoutput mapping is truly being learnt.
The second set of simulation experiments are with real
Electroencephalogram data collected at the Freiburg university
[26]. There are 6 channels in all the recordings because of use
of 6 electrodes. Each recording has been examined by an
expert encephalographer and occurrences of seizure identified.
Each recording covers a nearly continuous length of time that
included several Ictal as well as inter-Ictal intervals. To
shorten the run-times, middle of the inter-seizure intervals
were shortened (e.g. to a unit variance Gaussian random
number with a mean 262144 samples rather than original
intervals exceeding 921600 samples) but intervals leading
upto and the seizure intervals were included unmodified. For
purposes of the simulation experiment, for each recording,
measurement updates were performed through first several
seizures (“training” seizures), measurement updates were
stopped sometime before the occurrence of last few seizures
(“test” seizures) and seizure prediction was compared to
determination made by the expert encephalographer. Table
presents summary of these experiments.
These set of experiments were repeated for prediction ahead
into future by different amounts of time. The three sets
presented in the table corresponded to prediction ahead by 1,
1025 and 2049 instants. High ADRs were achieved for each of
these cases. To study whether, in spite of quantitatively similar
performance for prediction ahead into future by different
amounts of time, whether performance differed qualitatively.
Figure 16, Figure 17 and Figure 18 show close-up of when in
network determination recordings transition from no seizure to
seizure when predicting ahead by 1, 1025 and 2049 instants.

Figure 15: Evolution of weights

TABLE 1: SIMULATION EXPERIMENTS WITH EPILEPTIC SEIZURE DATA
EXPERIMENT
RECORD
RECORD
PREDICT
TRAIN
TEST
NAME
LENGTH
AHEAD
SAMPLES
SAMPLES
1
Pat003
2361858
1
2128306
233552
2
Pat003
2361858
1025
2128306
233552
3
Pat003
2302384
2049
2128304
174080

TABLE 2
REAL EPILEPTIC SEIZURE DATA TRAINING SEGMENTS
EXPERIMENT
TRAIN
TRAIN
SEIZURES
SAMPLES
MARKED
SEIZURE
1
3
111560
2
3
111560
3
3
111560

EXPERIMENT

1
2
3

TABLE 3
REAL EPILEPTIC SEIZURE DATA TEST RESULTS
TRUE
SPC
TRUE
SEN
POSITIVE,
NEGATIVE,
FALSE
FALSE
NEGATIVE
POSITIVE
6224, 175
0.97
227010,
1
143
6308, 91
0.97
227106,
1
47
6224175
0.97
167538,
1
143

ADR

0.986011
0.992786
0.99

The blue curve corresponds to expert determination and the
green curve to network determination. Similarly Figure 19,
Figure 20 and Figure 21 show close-up of when in network
determination recordings transition from seizure to no seizure
when predicting ahead by 1, 1025 and 2049 instants.

Figure 14: Predicted versus noise corrupted actual data for the Mackey Glass
series

15

Figure 16: 1 instant ahead prediction of seizure onset – predicted (green)
versus expert determination (blue)

Figure 18: 2048 instants ahead prediction of seizure onset – predicted
(dots) versus expert determination (solid line)

Figure 19: 1 instant ahead prediction of seizure subsiding
Figure 17: 1024 instants ahead prediction of seizure onset – predicted
(green) versus expert determination (blue)

Figure 20: 1024 instants ahead prediction of seizure subsiding

16
APPENDIX
Output of visible node
⎛
⎞
y v p = f ⎜⎜∑ w v p uq y˜ uq ⎟⎟
⎝ q
⎠
where
is set to 1,

€

(33)
is the

th input and

the weight from input to the visible node
∂y v p
= f ʹ′' ( a) a=a w v p uq
vp
∂y uq

is

. Hence
(34)

and

∂y v
€

Figure 21: 2048 instants ahead prediction of seizure subsiding

€
VI. CONCLUSION
We presented the time-series prediction problem€ and
critically discussed several aspects. Some of the aspects
pertain to applications with stringent requirements such
€ as
Epileptic Seizure onset prediction. We discussed suitability of
Recurrent Neural Networks for time series prediction. We
discuss an architecture for Recurrent Neural Networks that
grows along three dimensions – number of units in layers,
number of hidden layers and extent in number of instants of
back-propagation. We highlighted the similarity between RNN
with greater than 1 instant back-propagation extent and higher
order Markov processes. We prove analytically the fact that
additional hidden layers improve the approximation error rate.
Application of theory of back-propagation gradient descent
leads to a training method called Gradient Descent Backpropagation Through Time-and-Space. We derive a Dynamic
Programming (DP) procedure employing matrix operations for
the training. DP and use of matrix operations makes the
procedure efficient particularly when using data-parallel
libraries and on data-parallel or scalar architectures. The
simulations studies present the geometry of the parameter
space and verify using synthetic data that the network learns
the temporal structure in the data so that parameters converge
while model output displays same dynamic behavior as the
system. Simulation further showed the method to attain nearperfect Average Detection Ratio on real epileptic seizure data.
We have several related ongoing work that are
enhancements to reduce cost of software and Integrated
Circuit implementations while preserving the very good
performance on the discussed as well as other applications.
The enhancements fit number of hidden layers, number of
€ to
units in hidden layers and extent of back-propagation
problem and consider local minima avoidance in gradient
descent training. Applications include Natural Language
Processing, Automatic Speech Recognition and music
analysis.

p

∂w v p uq

= f 'ʹ′ ( a) a= a y uq .

(35)

vp

Alternatively, above can be written in matrix notation as
,
(36)
∂y v
= f ʹ′( a) a=a Wvu and
(37)
v
∂y u
∂y v
= f ʹ′( a) a=a .
(38)
v
∂Wvu
Now

is the vector of inputs to the network,

vector

is the

with a 1 as the additional and the first element

(corresponding to the bias input),

is the matrix of weights

to visible layer units from the inputs,

f ʹ′( a) a=a

(diagonal) matrix of derivatives of vector
evaluated at

and

is the
v

wrt vector

is the matrix of partial

derivatives of outputs of visible€ layer units wrt inputs.
, derivative of vector
wrt matrix
is also a
matrix (and not a 3-dimensional array) as derivative of the
output of a unit in a layer wrt weights to other units in the
same layer are 0.
Note that representing the set of derivatives as matrices is
more than just a matter of compactness. Representing as
matrices eases implementation on scalar processors using
mathematical library functions by reducing the overheads
associated with library function calls and on vector processors
such as GPUs enables gains from the data parallelism.
However in the rest of the following derivation, for ease of
understanding, we continue to use scalar notation and employ
sum over elements of vector as necessary.
Output of hidden layer
node

y h NL , j

NV
⎛ N H
⎞
= f ⎜⎜∑ w h h − y h − + ∑ w hNL , j v p y v p ⎟⎟
N , j N ,p
NL , p
p=1
⎝ p=1 L L
⎠

where

is the output of hidden layer
the weight from hidden layer

previous instant to hidden layer
instant,

is the output of the hidden layer

the previous instant,

node
node

node

(39)
,

at the

at present
node

at

the weight from the hidden layer

17

N L node j to visible node p and y v p the output of visible

propagating T instants, the derivative of final output y z i wrt

node p .
Hence

the weight associated with input q at instant −T to input layer
unit p is given by
NH
∂y hN , p €
€∂y zi
L
= f' ( a) a=a ∑ w zi hN , p
(49)
L
z i€
∂w ( − T ) ( − T )
∂w
€
p=1
( − T) ( −T)

∂y
€ h NL , j

€

€

= f 'ʹ′ ( a) a= a

∂y h −

€

NL , p

€

= f 'ʹ′ ( a) a= a

h NL , j h N−L , p

∂y hN

w€h

−
NL , j h NL , p

,

(40)

€

∂y hNL , j
∂w

h NL , j

h NL , p

and

yh −

(41)

NL , p

€

L,j

= f ʹ′' ( a) a=a
yvp .
(42)
hNL , j
∂w hN , j v p
L
From (33) and (39),
⎛ N H
NV
⎛ N I
⎞⎞
y hN , j = f ⎜∑ w h h ( −1) y h ( −1) + ∑ w hN , j v p f ⎜⎜∑ w v p uq y uq ⎟⎟⎟ (43)
L
⎜
L
⎟
NL , j NL , p
NL , p
p=1
⎝ q=1
⎠⎠
⎝ p=1
Hence
∂y hN , j
L
= f ʹ′( a) a=a
w hN , j v p f ʹ′( a) a=a y uq
L
hN L , j
vp
∂w v p uq

€
€

€

= f ʹ′' ( a) a=a

€

∂y hN

L,j

∂y uq
€

hNL , j

∂y v p

w hN , j v p
L

hNL , j

∑w
p=1

NV

= f ʹ′( a) a=a

hNL , j

∑w
p=1

v

v f ʹ′( a) a=ah

hNL , j p

NL , p

w v p uq

∂y v p

hNL , j p

(44)

∂y uq

Output of hidden layer N L −1,
⎛ N H
⎞
⎟⎟
y hN −1, j = f ⎜⎜ ∑ w h
y
+
χ
hN− −1, p hN− L −1, p
L
⎝ p=1 NL −1, j L
⎠
€
where

€
€

uq

vp

layer unit outputs and unit input weights are assigned indices
so as to facilitate computation of the derivatives.
are the current final outputs,
hidden

are the
layer

€

χ = ∑ wh

N L −1, j

p=1

those with hidden layers
(45)
(51)

NV
⎛ N I
⎞
+ ∑ w hN , p vq f ⎜⎜∑ w vqur y ur ⎟⎟
L
N L ,q
q=1
⎝ r=1
⎠

NH

q=1

€

Hence
∂y hN −1, j

€

∂w v p ur

and finally those with the input layer

hNL , p f (ξ )

ξ = ∑ wh

−
N L , p hN L ,q

y h−

(46)

.

L

= f ʹ′( a) a=a

hNL −1, j

NH

€

∑w

.

p=1

€

∂y hN

L −1,j

∂w v p ur

€

hNL −1, j hNL , p

f ʹ′( a) a=a

hNL , p

w hN

∂y vq
v
L,p q

∑w

a=ah

N L −1, j p=1

hNL −1, j hNL , p

(47)

∂w v p ur

∂y hN

NH

= f ʹ′( a)

L,p

∂w v p ur

(48)

(23) and (26) can in fact be generalized so we may compute
the derivative of a unit output or a weight wrt to a unit output
(say y x1 ) at an arbitrary distance away in terms of the

€

derivative of the same weight or output wrt to (an) output(s)
y x2 distance 1 away, the weight connecting y x2 to y x1 and

€

derivative of the original unit output or weight wrt y x2 . For
example in a network that has N L hidden layers, we are back

€

€

€

€
€

€

unit
outputs,
are the current

visible layer unit outputs. Current outputs are followed by
instant
, instant
, …, instant
outputs where
is
the extent of back-propagation. Weights follow the unit
outputs. Current weights associated with output layer units
occupy indices
,
(50)

NH

€

uq

Above suggests a dynamic programming method for
computing the derivatives comprising of traversing the graph
from inputs to outputs, computing derivative of unit outputs
wrt distance 1 away unit outputs, distance 1 away weights and
distance 2 away weights. Derivatives wrt unit outputs at
distance greater than 1 and weights at distance greater than 2
are computed by taking weighted sums over corresponding
derivatives of distance 1 away unit outputs. The distances here
are partially over space and partially over time.
In the joint-state-and-parameter-vector
, the individual

current

∂w v p uq
NV

= f ʹ′( a) a=a

vp

(52)

As with outputs, current weights are followed by instant
through instant
weights.
Before computing new current outputs, the current through
instant
outputs are copied to instant
through instant
outputs. Similarly before updating, the current through
instant
weights are copied to instant
through
instant
weights.
However significant time-efficiency is achieved by noting
that no actual copying needs to be performed. The unit outputs
and weights sections are treated as
(instants) long
circular buffers and the current instant pointer is simply
advanced one instant!
The row and column indices of matrices of derivatives of
unit outputs and weights wrt upstream unit outputs and
weights are also similarly incremented circularly when
advancing to next instant. This makes it unnecessary to

18
physically copy weights. As operations of matrix-vector
product and vector-vector inner product involves summing
over products of pairs of elements and the sum operation is
order-independent, the matrix operations are not affected. This
makes the time-complexity when using back-propagation
extent greater than 1 same as the time complexity when using
back-propagation extent equal to 1.
With N O long network outputs vector y z and the
derivatives of weights wrt network outputs arranged in a
N O × N W matrix J , the updated weights are given
by w ← w − ηJ( yˆ z − y z ) .
€
€

€

€

Thus, the space complexity is
. With
€
a constant independent of network dimensions, this is
linear in network dimension. With access to previous outputs
and weights requiring
time, the time-complexity of
computations at each instant is

under normally

valid assumptions of

.

and

ACKNOWLEDGMENT
Piyush was a summer intern at the Nano-Electronics Research
Laboratory in the Electrical Engineering Dept of Purdue
University for a part of this work and is grateful for support
and encouragement received from Prof K. Roy and Dr S.
Raghunathan.
REFERENCES
[1]

I. Arasaratnam et. al., “Discrete Time Nonlinear Filering Using GaussHermite Quadrature,” Proc. of IEEE, vol. x, no. x, pp. 953-977, May,
2007.
[2] Andrew S. Barron, “Universal Approximation Bounds for
Superpositions of a Sigmoidal Function,” IEEE Transactions on
Information Theory, vol. 39(3), pp. 930-945, 1993.
[3] Belhouari S. B. et. al., “Gaussian Process for Nonstationary Time Series
Prediction,” in Computational Statistics and Data Analysis, xth ed. City
of Publisher, Country if not USA: Elsevier, 2004, ch. x, sec. x, pp. xxx–
xxx.
[4] Bengio Y. et. al., “Deep Learning of Representations for Unsupervised
and Transfer Learning,” Journal of Machine Learning Research (Proc.
Unsupervised and Transfer Learning Challenge and Workshop), Vol.
27, 2012.
[5] R. van der Merwe and E. Wan, “Sigma-point Kalman filters for
probabilistic inference in dynamic state-space models,” in Proceedings
of the Workshop on Advances in Machine Learning, City of Conf.,
Abbrev. State (if given), 2003, pp. xxxxxx.
[6] R. van der Merwe, “Sigma-point Kalman filters for probabilistic
inference in dynamic state-space models,” Ph.D. thesis, OGI School of
Science & Engineering, Oregon Health & Science University, City of
Univ., Portland, OR, USA, April 2004.
[7] Duchi et. al. , “Ergodic Subgradient Descent,” Allerton Conference on
Communication, Control and Computing, 2011.
[8] Gallant A. R. et. al., “There exists a neural network that does not make
avoidable mistakes,” Proceedings of the Second Annual IEEE
Conference on Neural Networks, Vol 2, I.657 – I.664, 1988.
[9] Greene B R et al, “A comparison of quantitative EEG features for neonatal seizure detection,” Clinical Neurophysiology, 2008.
[10] Hornik K. et. al., “Universal Approximation of an Unknown Function
and its Derivatives using a Multilayer Feedforward Networks,” Neural
Networks, Vol 3, 550-560, 1990.
[11] A. H. Jazwinski, “Stochastic Processes and Filtering Theory,”
Academic Press, 1970.
[12] S. J. Julier, J. K. Uhlmann, and H. F. Durrant-Whyte, “A new method
for the nonlinear transformation of means and covariances in filters and

[13]
[14]
[15]
[16]
[17]
[18]

[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]

estimators,” IEEE Transactions on Automatic Control, vol. 45(3), pp.
477–482, March 2000.
Kimball B. F., “Note on Computation of Orthogonal Predictors,” The
Annals of Mathematical Statistics, Vol. 24, No. 2, 1953.
F. Le Dimet et. al. “Second Order Information in Data Assimilation,”
Journal of the American Metrological Society, Vol. 130, No. 3, pp 629648, March 2002.
Loh P.-L. et. al. “High-dimensional regression with noisy and missing
data: Provable guarantees with non-convexity,” available at
http://arxive.org/abs/1104.4824, 2012.
Niedermeyer E et. al., “Electroencephalography: Basic Principles,
Clinical Applications, and Related Fields,” Liipicott Williams &
Wilkins, 2004.
Raghunathan S et al, “ A hardware-algorithm co-design approach to
optimize seizure detection algorithms for implantable applications,”
Elsevier Journal of Neuroscience Methods, 2010.
Rhew H-G et al, “A Wirelessly Powered Log-based Closed-Loop Deep
Brain Stimulation SOC with Two-Way Wireless Telemetry for
Treatment of Neurological Disorders,” Symposium on VLSI Circuits,
2012.
F. Sandblom et. al., “Marginalized Sigma Point Form,” , 2011.
Sha H et al, “A micro-controller-based implantable nerve stimulator
used for rats,” IEEE Engineering in Medicine and Biology Society
Conference, 2006.
Subasi A., “Classification of EEG signal using neural network and
logistic regression,” Computer Methods and Programs in Bio-Medicine,
Vol. 78, No. 2, 2005.
Sudderth E. B. et. al. “Nonparametric Belief Propagation,”
Communications of the ACM, Vol. 53, No. 10, pp 95-103, October
2010.
Sun F. T. et. al., “Responsive Cortical stimulation for the treatment of
epilepsy,” Neurotherapeutics, Vol. 5, 2008.
Werbos P. J., “Backpropagation Thru Time: what it does and how to do
it,” Proceedings of the IEEE, Vol.78, No. 10, 1990.
Yeager D et al. “A Fully-Integrated 10.5 uW Miniaturized (0.125mm2)
Wireless Neural Sensor,” Symposium on VLSI Circuits, 2012.
Seizure Prediction Project, http://epilepsy.uni-freiburg.de/friburgseizure-prediction-project.

Sharat C Prasad (M’87) received the B.Tech. degree in
Electronics and Communication Engineering from I.I.T.
Kaharagpur, India in 1983, M.E. in Computer Science and
Automation from I.I.Sc. Bangalore India in 1987 and Ph.D. in
Electrical Engineering from University of Texas at Dallas in
2005.
He was an Intern at the Tata Institute of Fundamental
Research during the summer of 1982, worked at the Indian
Subsidiary of International Computers Limited, U.K. from
1983 until 1985, at Texas Instruments from 1987 until 1997, at
Cisco Systems from 1998 until 2006, at Google Inc from 2007
until 2012 and is with Best3 Corp since 2012.
Dr Prasad is a member of IEEE and ACM. He has authored
or co-authored one book, over 5 journal, over 10 conference
and over 30 patents.

Piyush Prasad was born in Plano, Texas, in 1995. He is
expecting to receive his B.S. degree in Biomedical
Engineering from Washington University, St Louis, MO in
2017.
He was a summer intern at the Nano-electronics laboratory
of Dept of Electrical Engineering at Purdue University.
Mr Prasad was a recipient of several awards at a national
Robotics league and is the co-author of a patent.

19

