# Deep Unsupervised Learning Using Spike-Timing-Dependent Plasticity

Sen Lu, Abhronil Sengupta
School of Electrical Engineering and Computer Science
The Pennsylvania State University
University Park, PA 16802, USA
Email: {senlu<EMAIL_ADDRESS>
###### Abstract

Spike-Timing-Dependent Plasticity (STDP) is an unsupervised learning mechanism
for Spiking Neural Networks (SNNs) that has received significant attention
from the neuromorphic hardware community. However, scaling such local learning
techniques to deeper networks and large-scale tasks has remained elusive. In
this work, we investigate a Deep-STDP framework where a convolutional network
is trained in tandem with pseudo-labels generated by the STDP clustering
process on the network outputs. We achieve $24.56\%$ higher accuracy and
$3.5\times$ faster convergence speed at iso-accuracy on a 10-class subset of
the Tiny ImageNet dataset in contrast to a $k$-means clustering approach.

###### Index Terms:

Unsupervised Learning, Spiking Neural Networks, Spike-Timing-Dependent
Plasticity

##  I Introduction

With high-quality AI applications permeating our society and daily lives,
unsupervised learning is gaining increased attention as the cost of procuring
labeled data has been skyrocketing concurrently. The ever-more data-hungry
machine learning models usually require a humongous amount of labeled data,
sometimes requiring expert knowledge, to achieve state-of-the-art performance
today. Since manual annotation requires a huge investment of resources,
unsupervised learning is naturally emerging as the best alternative.

One of the most prominent unsupervised learning methods is clustering. The
main concept of clustering is to compress the input data (like images in the
case of computer vision problems) into lower dimensions such that the low-
dimensional features can be clustered into separable groups. The efficiency of
the sample clustering process improves with better representations of the
compressed features. Since the quality of features depends only on the
dimension reduction algorithm, the design and choice of the clustering method
are critical to the success of unsupervised learning. However, most real-world
tasks are not easily represented as separable low-dimensional points. Earlier
attempts include classical PCA reduction before clustering [1], while others
attempt to augment more features with “bags of features” [2]; but mostly
constrained to smaller tasks. Recent works like DeepCluster have explored
scaling of unsupervised learning approaches by incorporating the $k$-means
clustering algorithm with a standard Convolutional Neural Network (CNN)
architecture that can learn complex datasets such as ImageNet without any
labels [3]. Some works have also proven that pre-training the network, even
unsupervised, is beneficial to building the final model in terms of accuracy
and convergence speed [4, 5, 6].

The focus of this article, however, is on scaling unsupervised learning
approaches in a relatively nascent, bio-plausible category of neural
architectures - Spiking Neural Networks (SNNs). SNNs have been gaining
momentum for empowering the next generation of edge intelligence platforms due
to their significant power, energy, and latency advantages over conventional
machine learning models [7, 8]. One of the traditional mechanisms of training
SNNs is through Spike-Timing-Dependent Plasticity (STDP) where the model
weights are updated locally based on firing patterns of connecting neurons
inspired by biological measurements [9]. STDP based learning rules have been
lucrative for the neuromorphic hardware community where various emerging
nanoelectronic devices have been demonstrated to mimic STDP based learning
rules through their intrinsic physics, thereby leading to compact and
resource-efficient on-chip learning platforms [10]. Recent works have also
demonstrated that unsupervised STDP can serve as an energy-efficient hardware
alternative to conventional clustering algorithms [11].

However, scaling STDP trained SNNs to deeper networks and complex tasks has
remained a daunting task. Leveraging insights from hybrid approaches to
unsupervised deep learning like DeepCluster [3], we aim to address this
missing gap to enable deep unsupervised learning for SNNs. Further, while
techniques like DeepCluster have shown promise to enable unsupervised learning
at scale, the impact of the choice of the clustering method on the learning
capability and computational requirements remains unexplored. The main
contributions of the paper can therefore be summarized as follows:

(i) We propose a hybrid SNN-compatible unsupervised training approach for deep
convolutional networks and demonstrate its performance on complex recognition
tasks going beyond toy datasets like MNIST.

(ii) We demonstrate the efficacy of STDP enabled deep clustering of visual
features over state-of-the-art $k$-means clustering approach and provide
justification through empirical analysis by using statistical tools, namely
Fisher Information Matrix Trace, to prove that STDP learns faster and more
accurately.

(iii) We also provide preliminary computational cost estimate comparisons of
the STDP enabled Deep Clustering framework against conventional clustering
methods and demonstrate the potential of significant energy savings.

##  II Related Works

Deep Learning: Unsupervised learning of deep neural networks is a widely
studied area in the machine learning community [12, 13]. It can be roughly
categorized into two main methods, namely clustering and association. Among
many clustering algorithms, $k$-means [14], or any variant of it [15, 16], is
the most well-known and widely used method that groups features according to
its similarities. Its applications can be found in practice across different
domains [17, 18]. Other approaches focus on associations to learn data
representations which are described by a set of parameters using architectures
such as autoencoders [19, 20] (where the data distribution is learnt by
encoding features in latent space).

In more recent works, such unsupervised learning methods have been applied to
larger and more complex datasets [3], making them applicable to more difficult
problems. Further, recent advances in generative models have also provided
opportunities at mapping unlabeled data to its underlying distribution,
especially in the domain of image generation using Generative Adversarial
Network (GAN) [21] with reconstruction loss directly [22] or using the auto-
encoded latent space [23, 22, 24]. Dumoulin $et$ $al.$’s recent effort at
combining GAN and auto-encoder has demonstrated even better performance [22].

Bio-Plausible Learning: Visual pattern recognition is also of great interest
in the neuromorphic community [25, 26]. In addition to standard supervised
vision tasks, SNNs offer a unique solution to unsupervised learning - the STDP
learning method [9]. In this scheme, the neural weight updates depend only on
the temporal correlation between spikes without any guiding signals, which
makes it essentially unsupervised. While it offers a bio-plausible solution,
it is rarely used beyond MNIST-level tasks[9, 27, 28] and primarily used for
single-layered networks. Going beyond conventional STDP based learning, Lee et
al. [27] proposed an STDP-based pre-training scheme for deep networks that
greedily trained the convolutional layers’ weights, locally using STDP, one
layer at a time but limited only to MNIST. Similarly, in Ferre et al.’s work
[29], the convolutional layers were trained on CIFAR10 and STL-10 with
simplified STDP, but the layers were also trained individually with complex
mechanisms. Further, their works are also limited to shallow convolutional
architectures.

Our work explores a hybrid algorithm design based on a merger of the above two
approaches. Our proposed framework provides a global training signal for the
CNN using a straightforward and end-to-end STDP-based SNN implementation. We
demonstrate significant accuracy improvement and computation savings for
VGG-15 architecture on the Tiny ImageNet dataset in contrast to state-of-the-
art deep clustering approaches.

##  III Preliminaries

###  III-A Deep Clustering with k-means Algorithm

Deep Clustering [3] enabled unsupervised training of visual features primarily
relies on the ability of clustering algorithms like the $k$-means to group
together similar data points. $k$-means is a popular unsupervised algorithm
for separating data points into distinct clusters. Given a user-specified
value of $k$, the algorithm will find $k$ clusters such that each data point
is assigned to its nearest cluster. The vanilla implementation of the
$k$-means algorithm iteratively calculates the Euclidean distance between
points for comparison and updates the cluster centroids to fit the given
distribution.

Deep Clustering utilizes the traditional CNN architecture to obtain the
features to be used for clustering. The reason behind this feature reduction
choice hinges upon the fact that a randomly initialized and untrained CNN
outperforms a simple multilayer perceptron network by a considerable margin
[30]. Driven by this observation, the main idea behind this framework is to
bootstrap the better-than-chance signal to teach the network and learn the
features. This teaching signal is transformed into a ‘pseudo-label’ so that
the network can learn from it. The ‘pseudo-labels’ which may or may not be the
same as the ground truth labels reflect the direction that the network weights
should be updated. By doing so, the feature extraction layers may become
slightly better at recognizing certain features and thereby producing more
representative features. The improved features can ideally be more separable,
thereby generating higher quality ‘pseudo-labels’. By repeating this process
iteratively, the CNN should ideally converge by learning the ‘pseudo-labels’
[3].

Note that the CNN layers used for feature-reduction purposes can be converted
into SNN layers with various methods as shown in many recent studies [31, 7,
32, 33, 34], or trained from scratch using backpropagation through time (BPTT)
[35, 36] which opens up the potential for adopting the entire feature-
reduction in a low-power neuromorphic setting. In this work, we therefore do
not focus on the CNN-SNN conversion and train it by backpropagation without
unrolling through time.

###  III-B STDP Enabled Neuromorphic Clustering

STDP is an unsupervised learning mechanism that learns or unlearns neurons’
synaptic connections based on spike timings [37]. In particular, the synaptic
connection is strengthened when the post-synaptic neuron fires after the pre-
synaptic neuron, and the connection is weakened if the post-synaptic neuron
fires before the pre-synaptic neuron. The intuition behind STDP follows
Hebbian learning philosophy where neurons that are activated together and
sequentially are more spatio-temporally correlated and thus form a pattern,
and vice versa. This learning rule enables the encoding of complex input
distributions temporally without the need for guiding signals such as the
label. The weights of the neuronal synapses are updated based on spike timings
[9] as follows:

$\Delta w=\begin{cases}A_{+}e^{-\frac{\Delta t}{\beta_{+}}},&\text{if}\ \Delta
t>0\\\ -A_{-}e^{\frac{\Delta t}{\beta_{-}}},&\text{if}\ \Delta t<0\end{cases}$
(1)

where, $w$ is the weight, $A_{+/-}$ are the learning rates, $\Delta t$ is the
exact time difference between post-neuron and pre-neuron firing and
$\beta_{+/-}$ are the time-constants for the learning windows. In practical
implementations, the exact spike timing is usually replaced with a spike trace
(see Section IV-B) that decays over time to reduce memory storage for STDP
implementation [38].

STDP training is predominantly explored in Winner-Take-All networks in
literature which consists of an excitatory layer of neurons with recurrent
inhibitory connections [9] (see “STDP Enabled SNN for Clustering” sub-panel in
Fig. 2). Such connections create a mechanism called ‘lateral inhibition’ where
activated neurons inhibit other neurons’ activities and therefore assist the
activated neurons to accentuate the learning process of its weights. To
prevent any neuron from dominating the firing pattern, the second key
mechanism is ‘homeostasis’ which balances the overall activities of the
neurons. Homeostasis prevents neurons from runaway excitation or total
quiescence. One popular way to achieve this is through adaptive and decaying
thresholding in which after every firing event, the firing threshold increases
such that the firing neuron requires higher membrane potential to fire again
in the future. Consequently, this will provide opportunities for other neurons
in the network to fire and learn the synaptic weights. The critical balance of
these two mechanisms ensures stable learning of the SNN. Fig. 1 shows an
example of STDP-trained weights of the excitatory neuron layer of an SNN where
representative digit shapes are learnt without any label information for the
MNIST dataset [39]. Each neuron in the network represents a cluster. By
running inferences on the STDP network, we can cluster the inputs according to
their corresponding most activated neuron. The learnt weights of each neuron
is equivalent to the centroid of the cluster represented by that neuron.

Figure 1: STDP learns generic features of input patterns (MNIST dataset) in
the excitatory synapses of the Winner-Take-All network. Each neuron represents
a cluster and its learnt weights represent the corresponding cluster centroid.

##  IV Methods

###  IV-A Proposed Deep-STDP Framework

Figure 2: Overall structure of Deep-STDP: The ConvNet compresses the input
images to a lower-dimensional feature vector which is mapped by STDP
clustering to a pseudo-label. The ConvNet is subsequently trained through
backpropagation using the pseudo-labels.

As mentioned previously, the convolutional layers of the network compress the
input images to a lower dimensional feature space as a one-dimensional vector.
In abstract terms, the framework solves the following optimization problem
[3]:

$\displaystyle\min_{w\in\mathbb{R}^{d\times
k}}\frac{1}{N}\sum_{n=1}^{N}\min_{y_{n}\in{\\{0,1\\}}^{k}}||f_{\theta}(img_{n})-w_{y_{n}}||_{1}$
(2) $\displaystyle\text{ such that }y_{n}^{\intercal}1_{k}=1$

where, $N$ is the total number of training samples, $y_{n}$ is the $n$-th
optimal neuron assignment encoded as a one-hot vector, $f_{\theta}$ is the
ConvNet forward pass output parameterized by its weights $\theta$, $img_{n}$
is the $n$-th input sample, $w_{y_{n}}$ is the STDP-learnt synaptic weight map
of the most activated neuron, $d$ is the feature dimension of the ConvNet
output and $k$ is the number of neurons/clusters in the network. By minimizing
the difference between the weights of the neurons and the patterns of the
features, we can obtain an SNN that generates optimal assignments of $y_{n}$
parameterized by weights $w$, which act as the pseudo-labels for our
algorithm.

With the pseudo-labels, the network training can be accomplished through the
standard minimization problem of network loss which can be described by:

$\min_{\rho,\theta}\frac{1}{N}\sum_{n=1}^{N}\mathcal{L}(g_{\rho}(f_{\theta}(img_{n})),y^{*}_{n})$
(3)

where, $\theta,\rho$ are parameters of the ConvNet $f_{\theta}(\cdot)$ and
classifier $g_{\rho}(\cdot)$ respectively, $\mathcal{L(\cdot)}$ is the loss
function, $img_{n}$ again is the $n$-th image input, $y^{*}_{n}$ is the $n$-th
optimal pseudo-label for this iteration.

However, SNNs only accept discrete spikes as input and therefore the ConvNet
feature outputs in floating-point representation (after appropriate pre-
processing like PCA reduction and $l_{2}$-normalization [3]) are subsequently
rate encoded by a Poisson spike train generator, where the feature values are
used as the Poisson distribution rate and sampled from the respective
distribution. At the end of the pseudo-label assignment, the STDP enabled SNN
resets for the next iteration. This is intuitive since after the ConvNet
weight update process, the feature distribution gets shifted and hence a new
set of neuron/cluster weights should be learnt by the STDP framework.
Algorithms 1-2 describe the overall structure of the proposed Deep-STDP
framework shown in Fig. 2.

Algorithm 1 DeepSTDP

1:$k=100$ // # of expected clusters$\times$10

2:$epochs\geq 200$ // training epochs limit

3:Initialize $ConvNet$ randomly

4:for $1$ to $epochs$ do

5: $ConvNet\Rightarrow$ Remove last layer

6: Initialize $DatasetX$ // Training images

7: for $img\in DatasetX$ do

8: Run $ConvNet$ on $img$ $\Rightarrow features$ // Inference

9: end for

10: Rate encodes $features$$\Rightarrow featureSet$

11: STDP($featureSet$) $\Rightarrow labels_{pseudo}$ // Train the SNN and
generate pseudo-labels using Algo. 2

12: $ConvNet\Rightarrow$ Adds back last layer

13: Train $ConvNet$ on $DatasetX$ and $labels_{pseudo}$ // Standard learning
through backpropagation

14:end for

Algorithm 2 STDP

1:$k=100$ // # of neurons

2:Table I parameters // See Table I for details

3:$timesteps=100$ // Simulation steps

4:$featureSet$ // Features to be trained on

5:Initialize $w_{+},w_{-}$ randomly

6:Initialize $k$ neurons with resting potential $V_{rest}$

7:Initialize $w_{inh}$ to -1 // Fixed inhibitory weights

8:$l=0$ // Refractory counter initialized for all neurons

9:$\epsilon=0$ // No extra threshold initially

10:for $(feature\in featureSet)$, $(t\Rightarrow 1$ to $timesteps)$ do

11: Encode $feature\Rightarrow s^{pre}_{+},s^{pre}_{-}$

12: Resets $\tau^{pre}_{+},\tau^{pre}_{-}$ to $\tau_{o}$ based on spikes

13: $\tau^{pre}_{+}=$ $exp(\frac{1}{\tau_{decay}})\tau^{pre}_{+}$

14: $\tau^{pre}_{-}=$ $exp(\frac{1}{\tau_{decay}})\tau^{pre}_{-}$

15: $\tau^{post}=$ $exp(\frac{1}{\tau_{decay}})\tau^{post}$ // Decay spike
traces

16: $\epsilon=exp(\frac{1}{\epsilon_{decay}})\epsilon$ // Decay adaptive
threshold

17: Decay potential using Eq. 4 $\Rightarrow V^{t}_{exc}$

18: $l=max(l-1,0)$ // Count down refractory period

19: for Neurons with $(l==0)$ do // Non-refractory

20: $V^{t}_{exc}\leftarrow$Update using Eq. 5

21: Fire spikes when $(V^{t}_{exc}>V_{thr}+\epsilon)\Rightarrow s^{post}$

22: for Neurons with ($V^{t}_{exc}>V_{thr}+\epsilon$) do

23: $l=L$ // Reset refractory period

24: $\epsilon=\epsilon+\alpha$ // Add adaptive threshold $\alpha$

25: Resets $\tau^{post}$ to $\tau_{o}$ based on spikes

26: end for

27: end for

28: Update using Eq. 7 and 8 $\Rightarrow\Delta w_{+},\Delta w_{-}$

29: $V^{t}_{exc}$ excluding itself
$\leftarrow\text{$s^{post}$}\cdot\text{$w_{inh}$}$ // Lateral inhibition

30:end for

31:return $labels_{pseudo}\Rightarrow$ most activated excitatory neuron IDs

###  IV-B STDP Enabled SNN for Clustering

Clustering in the SNN is mediated through the temporal dynamics of Leaky-
Integrate-Fire neurons in the excitatory layer. In the absence of any spiking
inputs, the membrane potential of neurons in the excitatory layer is
represented by $V_{exc}$ at timestep $t$, or simply $V_{exc}^{t}$. It
initializes with $V_{exc}^{t=0}=V_{rest}$ and decays as,

$V_{exc}^{t}=V_{rest}+exp(\frac{1}{V_{decay}})(V_{exc}^{t-1}-V_{rest})$ (4)

where, $V_{rest}$ is the resting potential and $V_{decay}$ is the potential
decay constant.

Prior works [9] on using SNNs for clustering have mainly dealt with simple
datasets without negative-valued features. This is in compliance with the
nature of STDP learning for positive valued spikes. However, in our scenario,
we consider negative valued spiking inputs as well in order to rate encode the
negative features provided as output of the ConvNet. In order to enable STDP
learning for negative inputs, we decompose the weight map into positive and
negative components to learn positive and negative spike patterns
respectively. Therefore, in presence of spikes, the excitatory layer’s neuron
membrane potential dynamics is updated as,

$V_{exc}^{t}\leftarrow\text{$s^{pre}_{+}$}\cdot\text{$w_{+}$}+\text{$s^{pre}_{-}$}\cdot\text{$w_{-}$}$
(5)

where, the membrane potential is denoted by $V^{t}_{exc}$ at timestep $t$, and
the input spikes and pre-synaptic weights are represented by $s^{pre}$ and $w$
respectively (with their positive and negative counterparts). It is worth
mentioning here that pre-neurons refer to the input neurons and post-neurons
refer to the excitatory layer neurons since the synapses joining them are
learnt by STDP.

Further, there is a refractory period $L$ parameter for every neuron which
will only allow execution of Eq. 4 and 5 if the refractory counter, $l$,
equals ‘0’. A spike will be generated when the membrane potential at the
current timestep is greater than the membrane threshold:

$s=\begin{cases}1&\text{if }(V^{t}_{exc}>V_{thr}+\epsilon)\text{ and }(l=0)\\\
0&\text{otherwise}\end{cases}$ (6)

where, $V_{thr}$ is the membrane threshold to fire a spike, $\epsilon$ is the
adaptive threshold parameter, $l$ is the refractory period counter which is
reset to $L$ upon a firing event and decays by 1 otherwise (thereby preventing
neurons from firing for $L$ timesteps after a spike). $V^{t}_{exc}$ resets to
$V_{reset}$ after firing a spike. The adaptive threshold parameter acts as a
balancer to prevent any neuron from being over-active (homeostasis) and is
incremented by parameter $\alpha$ upon a firing event and otherwise decays
exponentially at every timestep similar to Eq. 4:
$exp(\frac{1}{\epsilon_{decay}})\epsilon$. Every spike generated by a post-
neuron triggers a membrane potential decrement by an amount $w_{inh}$ for all
the other neurons except itself.

In the context of our implementation, we used the spike trace $\tau$ to
represent the temporal distance between two spikes. The spike trace value
peaks at its firing to $\tau_{o}$ and exponentially decay as time lapses:
$exp(\frac{1}{\tau_{decay}})\tau$. The weight updates are similarly separated
into positive and negative parts.
Pre-synaptic update:

$\begin{array}[]{l@{}l}\Delta w_{+}&=-\eta^{pre}(s^{pre}_{+}*\tau^{post})\\\
\Delta w_{-}&=\eta^{pre}(s^{pre}_{-}*\tau^{post})\end{array}$ (7)

Post-synaptic update:

$\begin{array}[]{l@{}l}\Delta w_{+}&=\eta^{post}(\tau^{pre}_{+}*s^{post})\\\
\Delta w_{-}&=\eta^{post}(\tau^{pre}_{-}*s^{post})\end{array}$ (8)

where, $\Delta w$ are the weight updates, $\eta^{pre},\eta^{post}$ are the
learning rates for pre- and post-synaptic updates respectively, $\tau$ is the
spike trace, and $s$ is the spiking pattern. Superscript $(^{pre})$,
$(^{post})$ indicates whether the trace or spike is from pre- or post-synaptic
neuron respectively, and the subscript $(_{+})$, $(_{-})$ indicates whether
the operation is for positive or negative input spikes. Note that the negative
$s^{pre}_{-}$ can be flipped easily by the distributive property of matrix
multiplication.

##  V Experiments and Results

###  V-A Datasets and Implementation

The proposed method was evaluated on the Tiny ImageNet dataset, which is a
center-cropped subset of the large-scale ImageNet dataset [40]. Unlike the
ImageNet 2012 dataset, which contains 1000 object categories, the Tiny
ImageNet dataset comprises of only 200 categories. Due to computation
constraints, we selected the first 10 classes from the Tiny ImageNet dataset
by the naming order and considered both the training and testing sets for
those corresponding classes in this work. All images were normalized to zero
mean and unit variance and shuffled to avoid any bias. We chose VGG15 as the
baseline network architecture with randomly initialized weights. Simulations
were conducted using the PyTorch machine learning library and a modified
version of the BindsNet toolbox [38] as the base platform for the experiments.
The results reported for the DeepCluster framework [3] were obtained without
any modification to the open-source codebase associated with the work, and its
hyperparameters were unchanged unless mentioned in this work. The ConvNet
learning rate was set to $1e-2$ and the number of clusters was set to 10 times
the number of classes (recommended as optimal in Ref. [3] and also found
optimal in the Deep-STDP framework). The training was performed for 200
epochs. All results obtained were run on 2 GTX 2080Ti GPUs and the associated
hyper-parameters used for the Deep-STDP framework can be found in Table I.

Parameter | Value | Description
---|---|---
$k$ | 100 | Number of neurons in the network
$V_{rest}$ | -65 | Resting membrane potential
$V_{reset}$ | -60 | Resetting membrane potential
$V_{decay}$ | 20 | Membrane potential decay constant
$V_{thr}$ | -52 | Threshold potential to fire a spike
$L$ | 5 | Refractory period
$\tau_{o}$ | 1.0 | Trace value when a spike fires
$\tau_{decay}$ | 100 | Trace decay constant
$\alpha$ | 0.45 | Adaptive threshold increment
$\epsilon_{decay}$ | $1e7$ | Adaptive threshold decay constant
$\eta_{pre}$ | $1e-3$ | Pre-synaptic learning rate
$\eta_{post}$ | $1e-6$ | Post-synaptic learning rate
$w_{inh}$ | $-1$ | Fixed inhibitory recurrent connection weights
$T$ | 400 | SNN simulation duration
TABLE I: Hyper-parameters for STDP training

Numerous cluster re-assignment frequencies were explored and ‘1’ (‘2’) was
found to be the optimal for Deep-STDP (DeepCluster), i.e. the pseudo-labels
were generated by passing the entire dataset once (twice) every epoch. Note
that this frequency represents the number of dataset iterations per epoch.
Following the evaluation method proposed by Zhang et. al [41], we froze all
network parameters and trained a linear layer at the output to evaluate the
efficiency of the model to capture the distribution of images in the training
set as well as its usage as a pre-trained model for general use cases. We
fixed the random seeds in each experiment such that the clustering process is
deterministic for a particular run. To prevent loss in generality, all
accuracy results reported here represent the average value over 5 independent
runs with different sets of random seeds.

###  V-B Evaluation Metrics

####  V-B1 Fisher Information

The Fisher information (FI) quantitatively measures the amount of information
retained in a statistical model after being trained on a given data
distribution [42]. Many prior works have used this metric to measure different
aspects of deep learning models including SNN models [43, 44]. Unlike prior
works, we use pseudo-labels to generate FI instead of ground-truth labels. FI
reflects the impact of weight changes on the ConvNet output. If the FI of
model parameters is small, we can conclude that the model’s learning
efficiency is poor since the weights can be pruned without affecting the
output, and vice versa. Therefore, this metric implicitly measures the quality
of the pseudo-labels.

Let us consider that the network tries to learn $y$ from a distribution $p$
parametrized by a set of weights $\theta$. Given samples $x$, the posterior
distribution is $p_{\theta}(y|x)$. The Fisher information matrix (FIM) is
defined as:

$F=\mathbb{E}_{x\sim X}\mathbb{E}_{y\sim p_{\theta}(y|x)}[\nabla_{\theta}\log
p_{\theta}(y|x)\nabla_{\theta}\log p_{\theta}(y|x)^{T}]$ (9)

where, $X$ is the empirical distribution of the actual dataset. However, the
exact FIM is usually too large to be computed directly and therefore the value
is usually approximated by its trace, which is given by:

$\mathrm{Tr}(F)=\mathbb{E}_{x\sim X}\mathbb{E}_{y\sim
p_{\theta}(y|x)}[||\nabla_{\theta}\log p_{\theta}(y|x)||^{2}]$ (10)

in which the expectations can be replaced by the averaged observation from the
dataset of $N$ samples:

$\mathrm{Tr}(F)=\frac{1}{N}\sum_{k=1}^{N}||\nabla_{\theta}\log
p_{\theta}(y|x)||^{2}_{2}$ (11)

where, $\mathrm{Tr}(F)$ is the trace of FIM, $\nabla$ is the partial
derivative operator. We follow the same implementation as the algorithm
specified in Ref. [44].

####  V-B2 Normalized Mutual Information

Further, following the Deep Clustering work [3], we also measured the
Normalized Mutual Information (NMI) metrics to evaluate mutual information
between two consecutive assignments of the STDP-enabled SNN, given by Eq. 12.

$\text{NMI}(y^{p},y^{p-1})=\frac{I(y^{p};y^{p-1})}{\sqrt{[H(y^{p})H(y^{p-1})]}}$
(12)

where, $y^{p},y^{p-1}$ are label assignments for epoch $p-1$ and $p$
respectively, $I(\cdot)$ is the mutual information function, and $H(\cdot)$ is
the entropy function. Since the assignments $y^{p},y^{p-1}$ are consecutive
and are generated from the same inputs, a high NMI value indicates a high
correlation between the two sets of assignments as well as stable assignments
of the pseudo-labels.

###  V-C Performance Evaluation

Figure 3: Average accuracy over 5 independent runs of Deep-STDP and
DeepClustering frameworks for the 10-class subset of Tiny ImageNet dataset.
Figure 4: Fisher information metric comparison: Deep-STDP retains more
information over epochs. Pale curves are the actual data while darker curves
are the smoothed plots.

Fig. 3 demonstrates that Deep-STDP based unsupervised feature learning
significantly outperforms DeepCluster approach based on $k$-means clustering.
The superior quality of pseudo-labels generated by Deep-STDP is also explained
empirically by the FIM trace variation over the learning epochs (see Fig. 4).
While both algorithms perform similarly during the initial stages, the
accuracy and FIM trace start improving significantly for the Deep-STDP
approach over subsequent epochs. Performance evaluation metrics (NMI, FIM and
Accuracy) for the two approaches at the end of the training process are
tabulated in Table II.

Metric | DeepCluster [3] | Deep-STDP
---|---|---
NMI | 0.5741 | 0.8597
FIM | 0.1890 | 0.5441
Acc | 0.3240 | 0.5696
TABLE II: Evaluation Metrics Comparison

In addition to training an additional linear layer for numerical performance
analysis, we also visualized the convolutional filter activations of the CNN
trained using our proposed framework. We can observe from Fig. 5 that the
network forms distinct filters specialized for completely different visual
patterns in different layers without using any ground truth label information.
On the other hand, similar visualization performed on the DeepCluster trained
network yielded similar simple patterns in the shallow layers without any
complex patterns represented in the deeper layers, further substantiating the
efficacy of the Deep-STDP approach.

Figure 5: Deep-STDP filter activations of Gaussian random noise from layers 6,
12, 24, and 36. We have used the unit-level visualization method proposed in
Ref. [45].

###  V-D Computational Cost Estimation

While a detailed system level hardware analysis for the two approaches is
outside the scope of this work, we provide motivation for neuromorphic deep
clustering by performing a comparative analysis of the computational cost of
the two approaches.

####  V-D1 Cost of k-means Clustering

To find the new centroid of a particular cluster, the algorithm calculates the
averaged center of all the data points assigned to that cluster using the
following equation:

$c_{j}=\frac{1}{|C_{j}|}\sum{x_{i}\in C_{j}}$ (13)

where, $c_{j}$ is the averaged coordinates of the $j$-th centroid, $|C_{j}|$
is the number of data points assigned to that corresponding cluster, and
$x_{i}$ is the $i$-th data point. Subsequently, the algorithm calculates the
Euclidean distance between every data point and every centroid and assigns
each data point to the cluster with the shortest distance to its centroid. The
goal is to solve the optimization problem:

$\operatorname*{argmin}_{C}\sum_{j=1}^{k}\sum_{i=1}^{|C_{j}|}||x_{i}-c_{j}||^{2}_{2}$
(14)

where, $\operatorname*{argmin}_{C}$ solves for the optimal centroids and $k$
is the total number of clusters. The above two calculations will be repeated
until convergence is achieved or until a maximum number of iterations is
reached. Hence, the number of mathematical operations can be summarized as
follows:

  * •

Clustering Step: Compute the distance $||x_{i}-c_{j}||^{2}_{2}$ from every
point to every centroid and assign to $k$ clusters

  * •

Update Step: Re-center the centroids in new clusters by averaging over
$|C_{j}|$ for all clusters

  * •

Repeat $it$ times

To calculate the distance of a point $x_{i}$ from $c_{j}$:

$||x_{i}-c_{j}||^{2}_{2}=\sqrt{\sum_{m=1}^{d=256}(x_{im}-c_{jm})^{2}}$ (15)

where, $d$ is the number of dimensions in the feature. Hence, the number of
multiplications (the number of squaring operations) in order to calculate the
Euclidean distance is:

$[k\cdot d]\cdot it\cdot N$ (16)

and the number of addition operations involved is:

$[k\cdot(2d-1)+d]\cdot it\cdot N$ (17)

where, $k$ is the number of clusters, $N$ is the number of training samples,
and $it$ is the number of maximum iterations in the $k$-means algorithm. In
Eq. 17, the $k\cdot(d-1)$ component arises from the summation of individual
distance along each dimension while another $k\cdot d$ component arises from
the subtraction operation for distance calculation along each dimension. The
last $d$ component arises from updating the new cluster coordinates (which in
the worst case will iterate through all data points, see Eq. 13). Given the
cost of float ADD operation is 0.9pJ and float MULT operation is 3.7pJ in 45nm
CMOS process [46], we estimated the total computational cost in the clustering
process for every training epoch to be 14.1mJ (considering $it=20$).
Considering 175 epochs of DeepCluster training to reach peak accuracy, the
total computational cost is 2467.5mJ.

####  V-D2 Cost of STDP Clustering

In the STDP based clustering approach, the computations can be summarized into
the following parts:

  * •

Feedforward Step: Integrate input Poisson spike train through the synapses
connecting input and excitatory layer

  * •

Learning Step: Updating the excitatory layer weights based on pre- and post-
synaptic spiking activities

  * •

Inhibition Step: Updating the neuron membrane potential based on lateral
inhibitory connections

  * •

Repeat $T$ times

Although multiplication symbols were used in Algo. 2, computation with spike
signals can always be reduced to summation operation since the spike magnitude
is always ‘0’ or ‘1’ [7]. Further, the addition operation is conditional upon
the receipt of spikes, thereby reducing the computation cost by a significant
margin for a highly sparse spike train. For instance, the average spiking
probability per neuron per timestep in the excitatory layer of the network is
only 0.19%. Hence, the total number of addition operations can be summarized
as:

$[p_{input}\cdot|w^{exc}|+(p_{input}+p_{exc})\cdot|w^{exc}|+p_{exc}\cdot|w^{inh}|]\\\
\cdot T\cdot N$ (18)

where, $p_{input},p_{exc}$ are the average (per neuron per timestep averaged
over the entire training process) spiking probability of the input and
excitatory neuronal layer respectively, $|w^{exc}|$ is the number of synaptic
connections between the input and excitatory layer, either $|w_{+}|$ or
$|w_{-}|$ since the input can be either positive or negative, $|w^{inh}|$ is
the total number of inhibitory connections in the network, $T$ is the number
of timesteps used for the STDP training process, and $N$ is the number of
training samples.

It is worth mentioning here that we primarily focus on the computationally
expensive portions of both algorithms for these calculations. In Eq. 18, the
$p_{input}\cdot|w^{exc}|$ component arises from the feedforward propagation of
input spikes, $(p_{input}+p_{exc})\cdot|w^{exc}|$ component arises from the
learning step and $p_{exc}\cdot|w^{inh}|$ arises from the inhibition step.
Therefore, the total computational cost for Deep-STDP per epoch is 55.34mJ and
considering 50 epochs of training (iso-accuracy comparison as shown in Fig.
3), the total energy consumption is estimated to be 2767.2mJ - comparable to
the DeepCluster framework.

####  V-D3 System Level Cost Comparison:

We note that the STDP based framework does not change the computational load
of the clustering framework significantly. However, the computational load at
the system level will be also dependent on the computational load for feature
extraction in the ConvNet. For instance, Ref. [3] mentions a third of the time
during a forward pass is attributed to the clustering algorithm while the
remaining is attributed to the deep ConvNet feature extraction. Therefore, we
expect the Deep-STDP based framework to be significantly more resource
efficient than the DeepCluster based approach due to $3.5\times$ reduction in
the number of training epochs - equivalently reducing the ConvNet feature
extraction computational cost.

##  VI Conclusions

In conclusion, we proposed an end-to-end hybrid unsupervised framework for
training deep CNNs that can be potentially implemented in a neuromorphic
setting. We demonstrated significant benefits in terms of accuracy and
computational cost by leveraging bio-plausible clustering techniques for deep
unsupervised learning of visual features and substantiated our claims by
empirical analysis through statistical tools like Fisher Information and
Normalized Mutual Information. Our work significantly outperforms prior
attempts at scaling bio-inspired learning rules like STDP to deeper networks
and complex datasets. Future work can focus on further scaling of the approach
and delving deeper into the mathematical underpinnings of the superior
performance of STDP as a deep clustering mechanism.

## Acknowledgments

This material is based upon work supported in part by the U.S. Department of
Energy, Office of Science, Office of Advanced Scientific Computing Research,
under Award Number #DE-SC0021562 and the National Science Foundation grant CCF
#1955815 and by Oracle Cloud credits and related resources provided by the
Oracle for Research program.

## References

  * [1] C. Ding and X. He, “K-means clustering via principal component analysis,” in _Proceedings of the twenty-first international conference on Machine learning_ , 2004, p. 29.
  * [2] G. Csurka, C. Dance, L. Fan, J. Willamowski, and C. Bray, “Visual categorization with bags of keypoints,” in _Workshop on statistical learning in computer vision, ECCV_ , vol. 1, no. 1-22. Prague, 2004, pp. 1–2.
  * [3] M. Caron, P. Bojanowski, A. Joulin, and M. Douze, “Deep clustering for unsupervised learning of visual features,” in _Proceedings of the European conference on computer vision (ECCV)_ , 2018, pp. 132–149.
  * [4] A. Radford, L. Metz, and S. Chintala, “Unsupervised representation learning with deep convolutional generative adversarial networks,” _arXiv preprint arXiv:1511.06434_ , 2015.
  * [5] A. v. d. Oord, Y. Li, and O. Vinyals, “Representation learning with contrastive predictive coding,” _arXiv preprint arXiv:1807.03748_ , 2018\.
  * [6] A. Radford, J. Wu, R. Child, D. Luan, D. Amodei, I. Sutskever _et al._ , “Language models are unsupervised multitask learners,” _OpenAI blog_ , vol. 1, no. 8, p. 9, 2019.
  * [7] A. Sengupta, Y. Ye, R. Wang, C. Liu, and K. Roy, “Going deeper in spiking neural networks: Vgg and residual architectures,” _Frontiers in neuroscience_ , vol. 13, p. 95, 2019.
  * [8] M. Davies, A. Wild, G. Orchard, Y. Sandamirskaya, G. A. F. Guerra, P. Joshi, P. Plank, and S. R. Risbud, “Advancing neuromorphic computing with loihi: A survey of results and outlook,” _Proceedings of the IEEE_ , vol. 109, no. 5, pp. 911–934, 2021.
  * [9] P. Diehl and M. Cook, “Unsupervised learning of digit recognition using spike-timing-dependent plasticity,” _Frontiers in Computational Neuroscience_ , vol. 9, p. 99, 2015.
  * [10] A. Saha, A. Islam, Z. Zhao, S. Deng, K. Ni, and A. Sengupta, “Intrinsic synaptic plasticity of ferroelectric field effect transistors for online learning,” _Applied Physics Letters_ , vol. 119, no. 13, 2021.
  * [11] E. P. Frady, G. Orchard, D. Florey, N. Imam, R. Liu, J. Mishra, J. Tse, A. Wild, F. T. Sommer, and M. Davies, “Neuromorphic nearest neighbor search using intel’s pohoiki springs,” in _Proceedings of the neuro-inspired computational elements workshop_ , 2020, pp. 1–10.
  * [12] Y. Bengio, A. C. Courville, and P. Vincent, “Unsupervised feature learning and deep learning: A review and new perspectives,” _CoRR, abs/1206.5538_ , vol. 1, no. 2665, p. 2012, 2012.
  * [13] H. U. Dike, Y. Zhou, K. K. Deveerasetty, and Q. Wu, “Unsupervised learning based on artificial neural network: A review,” in _2018 IEEE International Conference on Cyborg and Bionic Systems (CBS)_. IEEE, 2018, pp. 322–327.
  * [14] S. Lloyd, “Least squares quantization in pcm,” _IEEE transactions on information theory_ , vol. 28, no. 2, pp. 129–137, 1982.
  * [15] K. Krishna and M. N. Murty, “Genetic k-means algorithm,” _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ , vol. 29, no. 3, pp. 433–439, 1999.
  * [16] D. Arthur and S. Vassilvitskii, “K-means++ the advantages of careful seeding,” in _Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms_ , 2007, pp. 1027–1035.
  * [17] H. Ng, S. Ong, K. Foong, P.-S. Goh, and W. Nowinski, “Medical image segmentation using k-means clustering and improved watershed algorithm,” in _2006 IEEE southwest symposium on image analysis and interpretation_. IEEE, 2006, pp. 61–65.
  * [18] K.-j. Kim and H. Ahn, “A recommender system using ga k-means clustering in an online shopping market,” _Expert systems with applications_ , vol. 34, no. 2, pp. 1200–1209, 2008.
  * [19] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” _nature_ , vol. 323, no. 6088, pp. 533–536, 1986.
  * [20] G. E. Hinton and R. R. Salakhutdinov, “Reducing the dimensionality of data with neural networks,” _science_ , vol. 313, no. 5786, pp. 504–507, 2006\.
  * [21] R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer, “High-resolution image synthesis with latent diffusion models,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 10 684–10 695.
  * [22] P. Bojanowski, A. Joulin, D. Lopez-Paz, and A. Szlam, “Optimizing the latent space of generative networks,” _arXiv preprint arXiv:1707.05776_ , 2017.
  * [23] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” _arXiv preprint arXiv:1312.6114_ , 2013.
  * [24] J. Masci, U. Meier, D. Cireşan, and J. Schmidhuber, “Stacked convolutional auto-encoders for hierarchical feature extraction,” in _Artificial Neural Networks and Machine Learning–ICANN 2011: 21st International Conference on Artificial Neural Networks, Espoo, Finland, June 14-17, 2011, Proceedings, Part I 21_. Springer, 2011, pp. 52–59.
  * [25] P. U. Diehl, D. Neil, J. Binas, M. Cook, S.-C. Liu, and M. Pfeiffer, “Fast-classifying, high-accuracy spiking deep networks through weight and threshold balancing,” in _2015 International joint conference on neural networks (IJCNN)_. ieee, 2015, pp. 1–8.
  * [26] E. Neftci, S. Das, B. Pedroni, K. Kreutz-Delgado, and G. Cauwenberghs, “Event-driven contrastive divergence for spiking neuromorphic systems,” _Frontiers in neuroscience_ , vol. 7, p. 272, 2014.
  * [27] C. Lee, P. Panda, G. Srinivasan, and K. Roy, “Training deep spiking convolutional neural networks with stdp-based unsupervised pre-training followed by supervised fine-tuning,” _Frontiers in Neuroscience_ , vol. 12, 2018.
  * [28] D. Liu and S. Yue, “Event-driven continuous stdp learning with deep structure for visual pattern recognition,” _IEEE Transactions on Cybernetics_ , vol. 49, no. 4, pp. 1377–1390, 2019.
  * [29] P. Ferré, F. Mamalet, and S. J. Thorpe, “Unsupervised feature learning with winner-takes-all based stdp,” _Frontiers in computational neuroscience_ , vol. 12, p. 24, 2018.
  * [30] M. Noroozi and P. Favaro, “Unsupervised learning of visual representations by solving jigsaw puzzles,” in _Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part VI_. Springer, 2016, pp. 69–84.
  * [31] R. Midya, Z. Wang, S. Asapu, S. Joshi, Y. Li, Y. Zhuo, W. Song, H. Jiang, N. Upadhay, M. Rao _et al._ , “Artificial neural network (ann) to spiking neural network (snn) converters based on diffusive memristors,” _Advanced Electronic Materials_ , vol. 5, no. 9, p. 1900060, 2019.
  * [32] S. Lu and A. Sengupta, “Exploring the connection between binary and spiking neural networks,” _Frontiers in neuroscience_ , vol. 14, 2020.
  * [33] ——, “Neuroevolution guided hybrid spiking neural network training,” _Frontiers in neuroscience_ , vol. 16, 2022.
  * [34] H. Gao, J. He, H. Wang, T. Wang, Z. Zhong, J. Yu, Y. Wang, M. Tian, and C. Shi, “High-accuracy deep ann-to-snn conversion using quantization-aware training framework and calcium-gated bipolar leaky integrate and fire neuron,” _Frontiers in Neuroscience_ , vol. 17, p. 1141701, 2023.
  * [35] G. Bellec, D. Salaj, A. Subramoney, R. Legenstein, and W. Maass, “Long short-term memory and learning-to-learn in networks of spiking neurons,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [36] N. Rathi and K. Roy, “DIET-SNN: Direct input encoding with leakage and threshold optimization in deep spiking neural networks,” _ArXiv_ , vol. abs/2008.03658, 2020.
  * [37] N. Caporale and Y. Dan, “Spike timing–dependent plasticity: a hebbian learning rule,” _Annu. Rev. Neurosci._ , vol. 31, pp. 25–46, 2008.
  * [38] H. Hazan, D. J. Saunders, H. Khan, D. Patel, D. T. Sanghavi, H. T. Siegelmann, and R. Kozma, “Bindsnet: A machine learning-oriented spiking neural networks library in python,” _Frontiers in Neuroinformatics_ , vol. 12, p. 89, 2018\.
  * [39] L. Deng, “The mnist database of handwritten digit images for machine learning research,” _IEEE Signal Processing Magazine_ , vol. 29, no. 6, pp. 141–142, 2012.
  * [40] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “ImageNet: A large-scale hierarchical image database,” in _Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on_. IEEE, 2009, pp. 248–255.
  * [41] R. Zhang, P. Isola, and A. A. Efros, “Colorful image colorization,” in _Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part III 14_. Springer, 2016, pp. 649–666.
  * [42] S.-i. Amari and H. Nagaoka, _Methods of information geometry_. American Mathematical Soc., 2000, vol. 191.
  * [43] R. Karakida, S. Akaho, and S.-i. Amari, “Universal statistics of fisher information in deep neural networks: Mean field approach,” in _The 22nd International Conference on Artificial Intelligence and Statistics_. PMLR, 2019, pp. 1032–1041.
  * [44] Y. Kim, Y. Li, H. Park, Y. Venkatesha, A. Hambitzer, and P. Panda, “Exploring temporal information dynamics in spiking neural networks,” _arXiv preprint arXiv:2211.14406_ , 2022.
  * [45] D. Erhan, Y. Bengio, A. Courville, and P. Vincent, “Visualizing higher-layer features of a deep network,” _University of Montreal_ , vol. 1341, no. 3, p. 1, 2009.
  * [46] S. Han, J. Pool, J. Tran, and W. Dally, “Learning both weights and connections for efficient neural network,” _Advances in neural information processing systems_ , vol. 28, 2015.


#  Chaotic attractor reconstruction using small reservoirs - the influence of
topology

Lina Jaurigue  Technische Universität Ilmenau<EMAIL_ADDRESS>
###### Abstract

Forecasting timeseries based upon measured data is needed in a wide range of
applications and has been the subject of extensive research. A particularly
challenging task is the forecasting of timeseries generated by chaotic
dynamics. In recent years reservoir computing has been shown to be an
effective method of forecasting chaotic dynamics and reconstructing chaotic
attractors from data. In this work strides are made toward smaller and lower
complexity reservoirs with the goal of improved hardware implementability and
more reliable production of adequate surrogate models. We show that a
reservoir of uncoupled nodes more reliably produces long term timeseries
predictions than complex reservoir topologies. We then link the improved
attractor reconstruction of the uncoupled reservoir with smaller spectral
radii of the resulting surrogate systems. These results indicate that, the
node degree plays an important role in determining whether the desired
dynamics will be stable in the autonomous surrogate system which is attained
via closed-loop operation of the trained reservoir. In terms of hardware
implementability, uncoupled nodes would allow for greater freedom in the
hardware architecture because no complex coupling setups are needed and
because, for uncoupled nodes, the system response is equivalent for space and
time multiplexing.

##  1 Introduction

Data driven models are needed in a wide range of applications and there are
many statistical and machine learning methods of developing them. Low training
costs and energy efficient approaches are paramount for wide spread use,
particularly in edge devices. For data driven timeseries forecasting,
particularly chaotic timeseries forecasting, reservoir computing has been
shown to perform well. Reservoir computing is a machine learning approach
wherein a dynamical system is driven with an input timeseries and the
responses of the system to the input are linearly combined to approximate the
desired output timeseries [1, 2]. In constrast to other machine learning
methods, the training is restricted to the linear output weights and can be
performed via linear regression. Due to the low training cost and the fact
that the dynamical system remains unchanged in the training procedure,
reservoir computing is well suited for hardware implementation.

In this work strides are made toward smaller and lower complexity reservoirs
which can lead to easier hardware implementation and interpretable surrogate
models. We show that a reservoir of 20 uncoupled nodes more reliably produces
better long term timeseries prediction performance than complex reservoir
topologies. These results have implications for the hardware implementability,
as uncoupled nodes are easier to construct and can equally be mutliplexed in
space and time, giving more freedom for possible hardware architectures.
Furthermore, to gain understanding about the connections between the reservoir
topology and both the short- and long-term prediction performance, error
measures which quantify the training accuracy, short-term predictions and the
accuracy of the attractor reconstruction are analysed, as well as properties
of the trained system. It is shown that an uncoupled reservoir topology leads
to a smaller spectral radius of the final trained system and that this
correlates with more stable attractor reconstruction.

The low training cost of reservoir computing has also lead to extensive
theoretical and numerical studies. These have typically used larger numbers of
reservoir nodes than the 20 used in this study. In [3] the closed loop
operation of a reservoir trained on one-step-ahead prediction for the Lorenz
and Kuramoto-Sivashinsky systems was investigated. The focus of this study was
the estimation of Lyaponov exponents of the Lorenz and Kuramoto-Sivashinsky
systems and reservoirs of 300 nodes were used. In [4] conditions for good
short and long term prediction were established. The reservoir used in this
study were sparse Erdős-Renyi networks of 2000 nodes. In [5] the authors used
a valid forecasting time measure, the correlation dimension and the largest
Lyaponov exponent to quantify short and long term prediction quality. They
investigate the influence of small world, Erdős-Renyi and scale free network
topologies, as well as the influence of the spectral radius of the original
system. Networks of 300 nodes were used. In [6] various reservoir topologies
were also compared for networks of 100 nodes, and Bayesian optimization was
applied for hyperparameter optimisation. Optimisation and validation
strategies for chaotic attractor reconstruction are investigated in [7]. In
[8] various block diagonal topologies were tested, with networks of 400-600
nodes. A block size of one is specifically not considered. The authors of [8]
compare the block diagonal topology with other reservoir computing related
approaches in [9] and find good performance compared with the SINDy approach
[10], nonlinear vector regression (referred to as ”next generation reservoir
computing” [11]) and reservoir computing using a Erdős-Renyi network.
Performance comparisons of reservoir computing with other machine learning and
statistical methods are given in [12]. Applications of trained chaotic
attractors using reservoir computing have been investigated. For example, in
[13] chaos synchronisation and cryptography applications are investigated and
in [14] it is shown that such systems can be used to infer unseen attractors.

Several hybrid reservoir computing approaches for chaotic timeseries
forecasting have been developed. One approach is to couple the reservoir at
the output layer with a so-called knowledge-based model [15]. This approach
requires an approximate model for the underlying dynamics of the system to be
reconstructed and has been shown to perform well, for example on atmospheric
modelling [16]. Often an knowledge based model is not available, therefore in
[17] a hybrid reservoir-SINDy approach was investigate. There the knowledge-
based model was replaced with a model produced from data using the SINDy
approach [10]. Another important topic is the influence of input noise on the
quality of the reconstruction of chaotic attractors. The authors of [18] have
shown that injection noise can lead to improved short- and long-term
predictions. And, in [19] the influence of input noise on the stability of
reconstructed attractors, as well as the frequency components whereof, are
investigated.

For many applications it is not only the reconstruction or forecasting of the
dynamics that are relevant, it is often also desirable to gain insights into
the underlying systems which are producing the data. This is difficult or even
impossible with most machine learning approaches as they do not produce
interpretable models. One key step toward an interpretable model is to greatly
restrict the size of the model. In terms of reservoir computing, this means
restricting the dynamical degrees of freedom of the system used as the
reservoir. A commonly used type of reservoir is a network of randomly
connected nodes with nonlinear activation functions. For such a reservoir, the
size of the reservoir is given by the number of nodes. In the context of
interpretable models, the dynamics of the final trained system are of
particular relevance, and these have so far not been thoroughly investigated.

In this work the influence of the reservoir topology on the ability to predict
the short and long term behaviour of the Lorenz chaotic attractor is
investigated, while restricting the number of reservoir nodes. The paper is
structured as follows. In Section 2 the reservoir computing appproach is
explain and the reservoir model and network topologies used in this work are
introduced. In Section 3 the timeseries prediction task is presented and open-
and closed-open predictions are explained. In Section 4 various error measures
are defined and simulation parameters are given. In Section 5 the results are
presented and discussed. Finally, the conclusions are given in Section 6.

##  2 Reservoir Computing

Figure 1: (a) Sketch of the reservoir computing approach. An input is fed into
the reservoir with fixed random weights ${\bf W}_{in}$ (grey arrows) and the
output is produced by linearly combining states of the reservoir with trained
output weights ${\bf W}$ (white arrows). (b) In open-loop operation externally
sourced input data is fed into the reservoir to produce the output sequence.
(c) In closed-loop operation the trained system is run autonomously by feeding
the output back in as the next input step.

In reservoir computing the response of a dynamical system to an input sequence
is sampled a number of times and then a linear weighted sum of these sampled
responses is used to make predictions or assign labels. In contrast to other
machine learning approaches, the training is restricted to the linear output
weights, as illustrated in Fig. 1a. Training the reservoir requires the
following steps:

  * •

Feed a sequence of inputs of length $K_{B}+K_{T}$ into the reservoir.

  * •

For each input step sample $N$ responses of the reservoir.

  * •

Disregard the responses to the first $K_{B}$ inputs, where $K_{B}$ is long
enough that the reservoir state no longer depends on the initial state.

  * •

Collect the $N$ responses to the remaining $K_{T}$ inputs into a time ordered
state matrix ${\bf S}\in\mathbb{R}^{K_{T}}\times\mathbb{R}^{N+1}$, adding a
bias column of ones. The entry $s_{k,j}$ of the state matrix then corresponds
to the $j^{\textrm{th}}$ sampled response of the reservoir to the
$K_{B}+k^{\textrm{th}}$ input.

  * •

Find the vector of weights ${\bf w}\in\mathbb{R}^{N+1}$ that minimises the
difference between the target output ${\bf y}\in\mathbb{R}^{K_{T}}$ and the
output ${\bf\hat{y}}={\bf S}{\bf w}\in\mathbb{R}^{K_{T}}$. The solution to
this problem is given by

${\bf w}=({\bf S}^{\textrm{T}}{\bf S}+\lambda{\bf I})^{-1}{\bf
S}^{\textrm{T}}{\bf y},$ (1)

using the Moore-Penrose pseudoinverse, where $\lambda$ is the Tikhonov
regularisation parameter and ${\bf
I}\in\mathbb{R}^{N+1}\times\mathbb{R}^{N+1}$ the identity matrix. For multiple
targets ${\bf Y}=\left[{\bf y_{1}},..,{\bf y_{q}}\right]$ the matrix of output
weights ${\bf W}=\left[{\bf w_{1}},..,{\bf
w_{q}}\right]\in\mathbb{R}^{N+1}\times\mathbb{R}^{q}$ must be calculated.

###  2.1 The reservoir

Many different dynamical systems are suitable reservoirs and have been
investigated both numerically and experimentally. These include optical [20,
21, 22, 23, 24], optoelectronic [25, 26], micromechanical [27, 28], quantum
reservoirs [29, 30, 31, 32]. In this study networks of the following form are
used as the reservoirs:

${\bf x}\left(k+1\right)=\tanh\left({\bf W}_{int}{\bf x}\left(k\right)+{\bf
W}_{in}{\bf i}_{in}\left(k+1\right)\right),$ (2)

where ${\bf x}\left(k\right)\in\mathbb{R}^{M}$ is the vector describing the
state of the $M$ network nodes at discrete time $k$, ${\bf
W}_{int}\in\mathbb{R}^{M}\times\mathbb{R}^{M}$ is the matrix of internal
reservoir node couplings (adjacency matrix), ${\bf
i}_{in}\left(k\right)\in\mathbb{R}^{p}$ is the vector of $p$ input signals at
time $k$ and ${\bf W}_{in}\in\mathbb{R}^{M}\times\mathbb{R}^{p}$ is the matrix
of input weights. The function $\tanh\left(\cdot\right)$ is the nonlinear
function which is applied element wise and commonly referred to as the
activation function. The $k^{\textrm{th}}$ row of the state matrix is filled
by sampling $N$ of the reservoir states ${\bf x}\left(k\right)$, with $N\leq
M$ and $k\in\left[0,K_{T}\right)$. For $N=M$ the elements of the state matrix
are $s_{k,j}=x_{j}\left(k\right)$ with $x_{j}\left(k\right)$
($j\in\left[0,..S\right)$) being the elements of the reservoir state vector
${\bf x}\left(k\right)$.

Figure 2: Reservoir network topologies corresponding to the coupling matrices
${\bf\tilde{W}}_{rand}$ (Random), ${\bf\tilde{W}}_{ring}$ (Ring) and
${\bf\tilde{W}}_{id}$ (Uncoupled).

Three reservoir coupling matrix types are compared in this study; a random
coupling matrix ${\bf\tilde{W}}_{rand}$ with entries chosen from a uniform
distribution on the interval $[0,1]$, an identity matrix ${\bf\tilde{W}}_{id}$
and a diagonal matrix ${\bf\tilde{W}}_{ring}$ with the entries shifted one
down from the main diagonal:

${\bf\tilde{W}}_{ring}=\begin{pmatrix}0&0&0&\cdots&0&1\\\ 1&0&0&\cdots&0&0\\\
0&1&0&\cdots&0&0\\\ \vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&0&0&\cdots&1&0\end{pmatrix}.$ (3)

In each case the matrices are rescaled such that ${\bf
W}_{int}=\rho_{R}{\bf\tilde{W}}a_{0}$, where $a_{0}$ is the scaling factor
needed such that $\rho_{R}$ is the spectral radius of ${\bf W}_{int}$. The
spectral radius is defined as

$\rho({\bf W})=\max\\{|\lambda_{1}|,...,|\lambda_{N}|\\},$ (4)

where $\lambda_{1},...,\lambda_{N}$ are the eigenvalues of the $N\times N$
matrix ${\bf W}$. For the ring and identity matrices $a_{0}=1$ and for the
random matrix $a_{0}$ depends on the randomly drawn values.

For ${\bf W}_{int}={\bf W}_{id}$ the nodes are uncoupled in the training
phase, i.e. each node only couples to its own state in the previous time step
and the node degree is zero. This coupling scheme is identical to delay-based
time-multiplexed reservoir computing when the feedback delay time of the
reservoir is equal to the input clocktime [33]. When ${\bf W}_{int}={\bf
W}_{ring}$, each node is coupled to one neighbour, thereby forming a
unidirectional ring. Sketches of the three coupling topologies are shown in
Fig. 2.

If the coupling matrices are filled with zeros, then the reservoir reduces to
an extreme learning machine (ELM), i.e. a single layer of uncoupled (also no
self-coupling) nodes with random input weights [34]. An example of a hardware
implemented ELM is presented in [35], where an array of microresonators is
used as a photonic ELM. In this limit, the system has no memory of past inputs
and the method can also be interpreted as nonlinear vector regression where
the feature vector contains various nonlinear transforms of the current input
step. In the results section this limit is included and corresponds to the
$\rho_{R}=0$ results (see Figs. 3-6 in Section 5).

##  3 Timeseries prediction

The aim of this study is to investigate the influence of the topology on the
ability of the reservoir to predict chaotic dynamics and act as a surrogate
model. To perform multi-step-ahead predictions and to act as a surrogate model
the reservoir is trained to predict a timeseries one step ahead and then in
the prediction phase the predicted output is used as the next input step. This
is the so-called closed-loop operation and is illustrated in Fig. 1c (open-
loop operation is shown in Fig. 1b). By coupling the reservoir output back in
as the next input, the system becomes autonomous. The desired result would
then be that the dynamics of the trained autonomous system emulate those of
the original system which produced the training data. In this study the Lorenz
system [36] is used to test to reservoir performance.

###  3.1 Lorenz 63 system

The Lorenz system [36] is given by

$\displaystyle\frac{dx}{dt}=c_{1}y-c_{1}x,\quad\frac{dy}{dt}=x(c_{2}-z)-y,\quad\textrm{and}\quad\frac{dz}{dt}=xy-
c_{3}z.$ (5)

With $c_{1}=10$, $c_{2}=28$ and $c_{3}=8/3$ this system exhibits chaotic
dynamics and has a Lyaponov exponent of $\lambda_{L}\approx 0.91$. We generate
the Lorenz timeseries using Runge-Kutta fourth order numerical integration
with a time step of $h=10^{-3}$. We use all three variables, sampled with a
step size of $\Delta t=0.1$ and rescaled to the range $\left[0,1\right]$, as
the reservoir input ${\bf i}_{in}(k)=\left[X(k),Y(k),Z(k)\right]$ ($X$, $Y$,
$Z$ indicate the rescaled versions of $x$, $y$ ,$z$). In the training phase
the target is given by ${\bf y}(k)=\left[X(k+1),Y(k+1),Z(k+1)\right]$. The
step size of $\Delta t=0.1$ is chosen, as it has been shown to be well suited
for the reconstruction of the Lorenz attractor [37, 38].

###  3.2 The trained autonomous system

After training, the reservoir output can be written in terms of the reservoir
state:

${\bf\hat{y}}\left(k\right)={\bf x}(k)^{\textrm{T}}{\bf W},$ (6)

where the reservoir state vector is extended by one term, $x_{N+1}=1$, to
include the bias. Labelling the entries of the output weight matrix (${\bf
W}$) which correspond to the bias terms as
$W_{N+1,i}=W{{}^{\textrm{bias}}_{i}}$, the reservoir outputs can be written as

$\hat{y}_{i}(k)=\sum_{j=1}^{N}x_{j}(k)W_{j}^{i}+W{{}^{\textrm{bias}}_{i}},$
(7)

for $i\in\left[1,q\right]$, where $q$ is the number of target sequences. Then
the next input is given by $I_{i}(k+1)=\hat{y}_{i}(k)$, i.e. ${\bf
i}_{in}(k+1)=\left[\hat{y}_{x}(k),\hat{y}_{y}(k),\hat{y}_{z}(k)\right]$

The reservoir then becomes an autonomous dynamical system described by

${\bf x}\left(k+1\right)=\tanh\left({\bf W}^{a}_{int}{\bf
x}\left(k\right)+{\bf b}\right),$ (8)

where ${\bf W}^{a}_{int}$ is the modified coupling matrix and ${\bf b}$ is the
vector of constant bias terms:

${\bf b}=\sum_{i=1}^{q}{\bf w}_{in}^{i}W_{\textrm{bias}}^{i},$ (9)

with ${\bf W}_{in}=[{\bf w}_{in}^{1},...,{\bf w}_{in}^{q}]$. The righthand
side of the final system now only dynamically depends on the previous states
of the nodes.

For the case of reconstructing a three dimensional dynamical system ($q=3$)
the new coupling matrix is given by

${\bf W}^{a}_{int}={\bf
W}_{int}+\begin{pmatrix}\sum_{i=1}^{3}W_{in,1}^{i}W_{1}^{i}&\sum_{i=1}^{3}W_{in,1}^{i}W_{2}^{i}&\cdots&\sum_{i=1}^{3}W_{in,1}^{i}W_{N}^{i}\\\
\sum_{i=1}^{3}W_{in,2}^{i}W_{1}^{i}&\sum_{i=1}^{3}W_{in,2}^{i}W_{2}^{i}&\cdots&\sum_{i=1}^{3}W_{in,2}^{i}W_{N}^{i}\\\
\vdots&\vdots&\vdots&\vdots\\\
\sum_{i=1}^{3}W_{in,N}^{i}W_{1}^{i}&\sum_{i=1}^{3}W_{in,N}^{i}W_{2}^{i}&\cdots&\sum_{i=1}^{3}W_{in,N}^{i}W_{N}^{i}\end{pmatrix},$
(10)

the entries of which depend on the initial reservoir, the input weights and
the trained output weights. We will relate spectral radius of this new
coupling matrix to the ability of the reservoir to operate as a surrogate
model. The spectral radius of ${\bf W}^{a}_{int}$ is calculated using Eq. (4)
and will be denoted $\rho_{a}$.

##  4 Methods

###  4.1 Performance measures

####  4.1.1 Normalised root mean squared error

The performance of the open-loop (see Fig. 1b) one-step-ahead prediction is
quantified using the normalised root mean squared error (NRMSE), defined as

$\textrm{NRMSE}=\sqrt{\frac{\sum_{k=1}^{K_{T}}\left(y_{k}-\hat{y}_{k}\right)^{2}}{K_{T}\textrm{var}\left({\bf
y}\right)}},$ (11)

where $y_{k}$ are the target values, $\hat{y}_{k}$ are the outputs produced by
the reservoir computer, $K_{T}$ is the number of testing steps (i.e. the
length of the vector ${\bf y}$) and $\textrm{var}\left({\bf y}\right)$ is the
variance of the target sequence.

####  4.1.2 Valid Prediction Time

We quantify the short-term closed-loop (see Fig. 1c) performance using a valid
prediction time (VPT) as in [17]. The valid prediction time is given by

$t_{\textrm{VPT}}=\textrm{max}\\{t|\delta_{\textrm{VPT}}<0.4\\},$ (12)

where $t$ is the time determined by the discretisation step $\Delta t$ of the
training data set and number of iterations from the beginning of the
prediction phase, and $\delta_{\textrm{VPT}}$ is a measure of the prediction
error defined as:

$\delta_{\textrm{VPT}}=\frac{|{\bf y}(k)-{\bf\hat{y}}(k)|^{2}}{\langle|{\bf
y}(k)-\langle{\bf y}(k)\rangle|^{2}\rangle}.$ (13)

Here $\langle\cdot\rangle$ indicates the time average over the length of the
target timeseries, ${\bf y}(k)=\left[X(k+1),Y(k+1),Z(k+1)\right]$ is the
vector of target values at discrete time $k$ and
${\bf\hat{y}}(k)=\left[\hat{X}(k+1),\hat{Y}(k+1),\hat{Z}(k+1)\right]$ is the
vector of predicted values.

Various other definitions of prediction times are used in the literature. For
example in [8] a forecast horizon is defined as the time at which the
difference in the predicted and target trajectories exceeds the standard
deviation of the target trajectory. Care must be taken when comparing these
measures as they can lead to very different times.

####  4.1.3 Attractor Deviation

To quantify how well the desired attractor is reconstructed we use an approach
similar to that used in [18]. We divided the 3-dimensional phase-space into
cubes of volume $dX\times dY\times dZ$. Each cube is assigned a value of one
or zero; zero if the trajectory never passes through the cube and one if the
trajectory passed through at least once. We do this for both the reconstructed
trajectory and the target trajectory. We then add together the absolute value
of the difference of the values assigned to each cube for the reconstructed
and target trajectories:

$\textrm{ADev}=\sum_{i}\sum_{j}\sum_{k}|(c^{true}_{ijk}-c^{pred}_{ijk})|,$
(14)

where $c^{true(pred)}_{ijk}$ is one if the trajectory passed through the
corresponding element of the phase space at least once and otherwise zero. The
ADev value depends on the size of the cubes, on the length of the timeseries
and on the target trajectory. We choose a cube dimension of $0.1\times
0.1\times 0.1$ and a timeseries length of 5000 time units ($5000\Delta t$),
which is about 450 Lyaponov times ($1/\lambda_{L}\approx 1.1$) for the Lorenz
63 system with the standard parameters used here.

####  4.1.4 Power Spectral Density

To compare the frequency components of the target and reconstructed
trajectories we calculate the power spectral density (PSD) of the $Z$
component. This is done using the timeseries over the prediciton horizon of
$5000\Delta t$, applying a Hamming windowing function and then calculating the
fast Fourier transform (FFT). The PSD is then given by

$S(f)=|\tilde{Z}(f)|^{2},$ (15)

where $\tilde{Z}(f)$ is the FFT of $Z(t)$ with $t=k\Delta t$. Finally, the
spectrum is smoothed by calculating a running average with a range of 20
points.

###  4.2 Simulation Parameters

For all tasks we rescale the input and target sequences to the range between
zero and one. Before beginning the training and testing phases, the reservoirs
are initialised by feeding in a sequence of 10000 inputs ${\bf i}_{in}(k)$. In
the training phase 10000 inputs are used while in the testing phase 5000
inputs are used. The lengths of the training and testing phases were not
optimised, as this was not the focus of this study.

In the next section results are presented in dependence of the spectral radius
of the initial reservoir $\rho_{R}$. The testing error measures are averaged
over 100 realisations of the random input weights ${\bf W}_{in}$ and coupling
weights ${\bf W}_{rand}$. For each $\rho_{R}$ value new random weights are
drawn and the coupling matrix is rescaled such that the spectral radius is
equal to the desired $\rho_{R}$. In each realisation the initial conditions
for the Lorenz target trajectory are chosen randomly for the prediction phase,
meaning that the start of the prediction phase is randomly distributed over
the Lorenz attractor.

##  5 Results and Discussion

Figure 3: NRMSE of the $X$, $Y$ and $Z$ components for one-step-ahead
prediction in open-loop operation as a function of spectral radius $\rho$ of
the coupling matrix. The mean over 100 realisations of the random weights and
Lorenz trajectories is depicted for the uncoupled (teal), random (orange) and
ring (grey) reservoirs. The error bars indicate the standard deviation.

To judge the forecasting capabilities of a reservoir, typically only the
performance in closed-loop operation is considered. We take a different
approach, and assess the quality of the one-step-ahead by calculating the
NRMSE on a test timeseries of 5000 steps in open-loop operation, i.e. the true
Lorenz timeseries is used as the input throughout the testing phase. We do
this to gain understanding about the influence of the training accuracy on the
closed-loop prediction performance. Figure. 3 shows the NRMSE for the $X$, $Y$
and $Z$ components of the Lorenz timeseries for the three reservoir
topologies. As a function of the spectral radius of the coupling martices,
$\rho_{R}$, the NRMSEs are very similar for all three topologies, with the
random network producing slightly lower errors. The lowest errors are achieved
near $\rho_{R}=0$, which is the ELM limit. For $\rho_{R}=0$ the reservoir has
no memory of past inputs, meaning that the predictions are based entirely on
the current input step. This lack of memory is not detrimental to the task
performance since the state of the Lorenz system at any point in time is fully
defined by the current $X$, $Y$ and $Z$ values and these are all being
provided as input. A study of the difference in reservoir performance for
partial and full input information in given in [39].

Figure 4: Percentage of the 100 realisations of the random weights and Lorenz
trajectories for which the resulting closed-loop (autonomous) trajectory; (a)
was bounded by $|X|,|Y|,|Z|<2$ for the entire prediction phase of 5000 steps,
(b) was bounded and remained oscillatory for the entire prediction phase of
5000 steps.

In the closed-loop (autonomous) operation, the predicted $X$, $Y$ and $Z$
values are fed back into the reservoir as the input. In this case three types
of resulting dynamics are found; diverging, oscillatory and fixed point. In
the diverging case the dynamic variables become much larger than the target
values and therefore the attractor is not successfully reconstructed. Since we
are interested in accurate timeseries prediction and attractor reconstruction,
it is not relevant in this study if the dynamics are truely diverging or if
the resulting system produces large oscillations. We therefore set an upper
limit of two and check if $|X|$, $|Y|$ and $|Z|$ remain below this value for
the entire prediction phase. Trajectories that remain below this limit we
refer to as bounded. In Fig. 4a the percentage of the realisations which were
bounded is depicted as a function of the spectral radius. Here, a clear
difference can be seen between the uncoupled reservoir and the random and ring
topologies. The uncoupled reservoir produces near 90% bounded trajectories
over a much large range of $\rho_{R}$. The bounded results can be separated
into those that remain oscillatory for the entire prediction phase and those
that collapse into a fixed point. The percentage of the results which remain
oscillatory for the complete prediction phase is shown in Fig. 4b. Here, a
clear difference is also evident between uncoupled reservoir and the random
and ring topologies. Figure 4 shows that the uncoupled reservoir produces
bounded results more reliably and that fewer of the results collapse into a
fixed point.

Figure 5: Closed-loop prediction performance: (a)-(c) The VPT for the three
reservoir topologies. (d)-(f) The attractor deviation ADev, with a reference
ADev calculated for two true Lorenz trajectories indicated by the dashed line.
(g)-(i) The time interval over which the trajectories were oscillatory. For
the ADev and run time values only the bounded trajectories were considered.
For all plots the median is plotted and the shaded regions show the
25${}^{\textrm{th}}$ to 75${}^{\textrm{th}}$ percentiles.

To evaluate the short term closed-loop prediction performance the valid
prediction time is plotted in Fig. 5a-c for the uncoupled, random and ring
topologies, respectively. For $\rho_{R}=0$ all three systems are identical,
and consequently show very similar performance as $\rho_{R}$ tends to zero.
The best performance in terms of the VPT is achieved near $\rho_{R}=0$ for all
three topologies. This correlates with the $\rho_{R}$ values for which the
lowest NRMSE errors are achieved in the open-loop operations, see Fig. 3. But,
in the case of the uncoupled reservoir, this does not correspond to the
$\rho_{R}$ values for which bounded and oscillatory trajectories are produced
most reliably, see Fig. 4. The finding that the highest VPT of to the Lorenz
task are achieved at $\rho_{R}=0$ are in agreement with the results presented
in [40]. In [40] reservoir hyperparameters are optimised for various tasks and
for the Lorenz task it is shown that the best prediction performance is
achieved when the adjacency matrix density is zero, which is equivalent to
$\rho_{R}=0$ in this study.

Figure 6: VPT for the uncoupled reservoir (orange) and for the true Lorenz
system (grey). For the true Lorenz VPTs the first prediction step of the
trained reservoir was used as the initial conditions the Lorenz equations and
these were then integrated to generate the predicted timeseries. The median is
plotted and the shaded regions show the 25${}^{\textrm{th}}$ to
75${}^{\textrm{th}}$ percentiles.

To judge the absolute performance it is useful to compare with other studies.
However, when making such comparisons differences, in the task implementation
and the error measures must be considered. In [17] the Lorenz task is
implemented with the same time discretisation as in this study and the same
definition for the valid prediction time is used. There VPTs of approximately
2 are achieved using a network of 100 nodes. By coupling the reservoir with a
SINDy model, this is improved to VPTs of about 3. In [15] a hybrid reservoir
of 50 nodes coupled with a knowledge based model, resulted in VPTs of
appproximately 1.5. The VPTs in [17] and [15] are larger than in this study,
however this was only achieved with larger node numbers and hybrid reservoir
methods.

As well as information about the performance relative to other reservoir
computing or machine learning approaches, we are also interested in the
performance relative to the Lorenz system itself. The largest Lyaponov
exponent of the Lorenz system indicates how quickly two trajectories on the
true attractor will diverge given a certain difference in the initial
conditions. Therefore, to obtain an objective benchmark with which the
absolute performance in terms of the VPT can be compared, we calculate the VPT
of the true Lorenz system given the error after the first prediction step (one
prediction step is made using the trained reservoir, then this prediction is
used as the initial conditions for the Lorenz system to generate a second
timeseries with which to compare the target trajectory). This comparison is
useful as it sets an upper limit on the VPT that can be expected given the
accuracy in the one-step-ahead training stage. This comparison is depicted in
Fig. 6 for the uncoupled reservoir. As is to be expected, the VPTs of the true
Lorenz system are slightly larger, however the 75${}^{\textrm{th}}$ percentile
of the reservoir results coincide with the median of the true Lorenz VPTs,
indicating considerable overlap in the performance.

Figure 7: Phase spaces plots (a,c,e) and the corresponding PSD plots (b,d,f)
for three example trajectories generated using the uncoupled reservoir in
closed-loop operation. The predicted results are shown in orange and the
target results in green. The trajectories in (c,d) and (e,f) correspond to
trajectories two and three in Fig. 8, respectively.

Figure 6 shows that even for the true Lorenz system, given differences in the
initial conditions of the order of the NRMSEs depicted in Fig. 3, trajectories
diverge after one to three Lyaponov times ($1/\lambda_{L}\approx 1.1$).
Therefore, to evaluate the quality of the attractor reconstruction over
hundreds of Lyaponov times we use the attractor deviation measure introduced
in Section 4.1.3. The attractor deviation (ADev) is a measure of the
difference in the 3-dimensional phase space which is traversed by the
trajectory of the target Lorenz timeseries and the predicted timeseries.
Figure 5d-f shows the attractor deviation for the uncoupled, random and ring
reservoirs. Only the realisations resulting in bounded trajectories were
considered. The horizontal dashed line is a reference attractor deviation of
two true Lorenz trajectories. The random and ring topologies show very similar
performance, and their large attractor deviations can be explained by the low
percentage of trajectories which do not collapse into a fixed point during the
prediction phase (as shown in Fig. 4b). The uncoupled reservoir, on the other
hand, shows quite low attractor deviation values corresponding to the
$\rho_{R}$ range for which a high percentage of the trajectories remain
oscillatory for the entire prediction phase.

Figure 8: Timeseries of the $X$, $Y$ and $Z$ components for two predicted
trajectories generated using the uncoupled reservoir in closed-loop operation.
The predicted results are shown in orange and the target results in green.
Trajectories two and three correspond to the trajectories in (c,d) and (e,f)
of Fig. 7, respectively.

Since the attractor deviation measure does not capture any temporal
information, it is also useful to compare the power spectral densities of the
predicted and true trajectories. Figure 7 shows phase-space plots, power
spectral densities and the attractor deviation for three example trajectories
generated using the uncoupled reservoir. Different random input weights and
initial conditions of the Lorenz system were used to generate the
trajectories. The trajectory in Fig. 7b,c collapses near the end of the
prediction phase and corresponds to trajectory two in Fig. 8. Figure 7e,f
corresponds to trajectory three in Fig. 8. As long as the trajectory remains
on the attractor, the spectral features are reproduced quite well and are
comparable with the performance achieved using a data-informed-reservoir
hybrid model in [17], even though a much smaller reservoir is used in this
study.

The Lorenz-like dynamics achieved with the trained reservoir appear to be
meta-stable. Meaning that the trajectories can appear Lorenz like for hundreds
of Lyaponov times and then suddenly collapse to a stable fixed-point, as
depicted in the timeseries shown in Fig. 8. Similiar behaviour has also been
reported in other works [5]. Comparing the results depicted in Fig. 3, Fig. 4
and Fig. 5 a trade-off between the short-term prediction performance and the
stable attractor reconstruction is evident.

Figure 9:  (a) Spectral radius of the trained autonomous system $\rho_{a}$ as
a function of the spectral radius of the original reservoir $\rho$. The solid
lines indicate the mean and the error bars give the standard deviation. (b)
Enlargement of the $\rho_{R}\in\left[0,0.5\right]$ range of (a).

For the reservoir to accurately reconstruct the Lorenz attractor, the Lorenz-
like dynamics of the trained system must be stable. For the uncoupled
reservoir a large percentage of the solutions appear to be stable or meta-
stable over the prediction range. However, for the random and ring topologies
most of the solutions diverge quickly, indicating that the desired solutions
either do not exist in the final autonomous system or are unstable. To work
towards understanding the difference in the results between the reservoir
topologies, the spectral radius of the trained autonomous systems $\rho_{a}$
is plotted in Fig. 9 as a function of the spectral radius of the original
reservoir. For all three reservoirs $\rho_{a}$ increases with $\rho_{R}$ and
for $\rho_{R}>0.5$ the large $\rho_{a}$ values correspond to short run times
(time before the system diverges or collapses to a fixed point), as shown in
Fig. 5g-i. For $\rho_{R}<0.35$, $\rho_{a}$ of the uncoupled system is smaller
than those of the random and ring reservoirs (see Fig. 9b for an enlargement
of the $\rho_{R}<0.5$ results). The region of low $\rho_{a}$ for the uncoupled
system corresponds to the regions of lower attractor deviation and longer run
times in Fig. 5d,g. These results indicate that the topology of the initial
reservoir plays a deciding role in determining the stability of the desired
solution in the final system. Specifically, less connectivity within the
network leads to improved results. These results are in agreement with the
findings of [8] where the authors showed improved performance for block
diagonal coupling matrices with low node degree (zero node degree was not
considered).

##  6 Conclusion

In this work the chaotic Lorenz attractor has been reconstructed using small
reservoirs of only 20 nodes. We have shown that the short term performance
depends on the error in the one-step-ahead prediction, but that the stable
reconstruction of the attractor is related to dynamical properties of the
trained autonomous system. By reducing the reservoir to a system of uncoupled
nodes, smaller spectral radii of the trained system can be achieved compared
with random or ring reservoirs. And, the parameter regions resulting in small
spectral radii for the trained uncoupled reservoir correspond to more reliable
reconstruction of the Lorenz attractor.

These results pave the way towards interpretable surrogate models. Due to the
small size of the reservoir and the uncoupled reservoir nodes it becomes
feasible to perform bifurcation analysis on the final system and additional
optimisation can be performed after the training phase without effecting all
reservoir nodes, for example pruning of selected nodes. Furthermore, hardware
implementability is improved, as complex network topologies are not required
and for a system of uncoupled nodes there is more freedom in the choice of
multiplexing.

The insights gained in this study into the relationships between training
accuracy, network topology and short- and long-term prediction performance can
also be applied to related machine learning approaches. For example, nonlinear
vector regression (referred to as ”next generation reservoir computing” [11])
has recently been shown to perform very well at chaotic timeseries
forecasting. Nonlinear vector regression could be thought of a network of node
degree zero (and no self-coupling) and would then be similar to the uncoupled
reservoir in this study.

In this study all three dynamical variables of the Lorenz system were provided
to the reservoir in the training stage. This meant that the reservoir did not
need to have memory for the one-step-ahead training to perform well. In
contrast, for tasks where only partial information is provided in the training
stage the reservoir must perform a delay embedding and therefore requires
memory [41, 42, 39]. Since the memory of the reservoir is determined by the
coupling topology and the spectral radius, different dependencies on these
properties can be excepted for partial information tasks compared with the
results for the full information Lorenz task presented in this study. In order
to restrict the reservoir topology to uncoupled nodes for partial information
tasks, but still incorporate the required memory, a possible solution could be
to apply memory augmentation methods at the levels of the input or output
layers. In [43] it was shown that adding task relevant timescales at the level
of the input or output layers can improve the reservoir performance for
partial information time series prediction tasks and decreases the sensitivity
to other reservoir hyperparameters.

## References

## References

  * [1] Jaeger H 2001 The ’echo state’ approach to analysing and training recurrent neural networks GMD Report 148 GMD - German National Research Institute for Computer Science
  * [2] Maass W, Natschläger T and Markram H 2002 Neural. Comput. 14 2531 ISSN 0899-7667
  * [3] Pathak J, Lu Z, Hunt B R, Girvan M and Ott E 2017 Chaos 27 121102 (Preprint https://doi.org/10.1063/1.5010300)
  * [4] Lu Z, Hunt B R and Ott E 2018 Chaos 28 061104 ISSN 1054-1500 (Preprint https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/1.5039508/14616258/061104_1_online.pdf) URL https://doi.org/10.1063/1.5039508
  * [5] Haluszczynski A and Räth C 2019 Chaos 29 103143 ISSN 1054-1500 (Preprint https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/1.5118725/14626476/103143_1_online.pdf) URL https://doi.org/10.1063/1.5118725
  * [6] Griffith A, Pomerance A and Gauthier D J 2019 Chaos 29 123108
  * [7] Racca A and Magri L 2021 Neural Netw. 142 252–268 ISSN 0893-6080
  * [8] Ma H, Prosperino D, Haluszczynski A and Räth C 2023 Chaos 33 ISSN 1054-1500 063130 (Preprint https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/5.0151290/17996145/063130_1_5.0151290.pdf) URL https://doi.org/10.1063/5.0151290
  * [9] Ma H, Prosperino D and Räth C 2023 Sci. Rep. 13 12970 ISSN 2045-2322 URL https://doi.org/10.1038/s41598-023-39886-w
  * [10] Brunton S L, Proctor J L and Kutz J N 2016 PNAS 113 3932–3937 (Preprint https://www.pnas.org/doi/pdf/10.1073/pnas.1517384113)
  * [11] Gauthier D J, Bollt E M, Griffith A and Barbosa W A S 2021 Nat. Commun. 12 5564 ISSN 2041-1723
  * [12] Gilpin W 2023 Phys. Rev. Res. 5(4) 043252 URL https://link.aps.org/doi/10.1103/PhysRevResearch.5.043252
  * [13] Antonik P, Gulina M, Pauwels J and Massar S 2018 Phys. Rev. E 98(1) 012215 URL https://link.aps.org/doi/10.1103/PhysRevE.98.012215
  * [14] Röhm A, Gauthier D J and Fischer I 2021 Chaos 31 103127 (Preprint https://doi.org/10.1063/5.0065813)
  * [15] Pathak J, Wikner A, Fussell R, Chandra S, Hunt B R, Girvan M and Ott E 2018 Chaos 28 041101 (Preprint 1803.04779)
  * [16] Arcomano T, Szunyogh I, Wikner A, Pathak J, Hunt B R and Ott E 2022 J. Adv. Model. Earth Syst. 14 e2021MS002712 e2021MS002712 2021MS002712 (Preprint https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2021MS002712) URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021MS002712
  * [17] Köster F, Patel D, Wikner A, Jaurigue L C and Lüdge K 2023 Chaos 33 073109 ISSN 1054-1500 URL https://doi.org/10.1063/5.0152311
  * [18] Zhai Z M, Kong L W and Lai Y C 2023 Phys. Rev. Res. 5(3) 033127 URL https://link.aps.org/doi/10.1103/PhysRevResearch.5.033127
  * [19] Wikner A, Harvey J, Girvan M, Hunt B R, Pomerance A, Antonsen T and Ott E 2024 Neural Networks 170 94–110 ISSN 0893-6080 URL https://www.sciencedirect.com/science/article/pii/S0893608023006159
  * [20] Brunner D, Soriano M C, Mirasso C R and Fischer I 2013 Nat. Commun. 4 1364
  * [21] Larger L, Baylón-Fuentes A, Martinenghi R, Udaltsov V S, Chembo Y K and Jacquot M 2017 Phys. Rev. X 7(1) 011015
  * [22] Röhm A, Jaurigue L C and Lüdge K 2019 IEEE J. Sel. Top. Quantum Electron. 26 7700108
  * [23] Sugano C, Kanno K and Uchida A 2020 IEEE J. Sel. Top. Quantum Electron. 26 1500409 ISSN 1558-4542
  * [24] Skontranis M, Sarantoglou G, Bogris A and Mesaritakis C 2022 Opt. Mater. Express 12 4047–4060 URL https://opg.optica.org/ome/abstract.cfm?URI=ome-12-10-4047
  * [25] Paquot Y, Duport F, Smerieri A, Dambre J, Schrauwen B, Haelterman M and Massar S 2012 Sci. Rep. 2 287
  * [26] Chen Y, Yi L, Ke J, Yang Z, Yang Y, Huang L, Zhuge Q and Hu W 2019 Opt. Express 27 27431
  * [27] Dion G, Mejaouri S and Sylvestre J 2018 J. Appl. Phys. 124 152132
  * [28] Yoshimura T, Haga T, Fujimura N, Kanda K and Kanno I 2023 Japanese Journal of Applied Physics 62 SM1013 URL https://dx.doi.org/10.35848/1347-4065/ace6ab
  * [29] Chen J, Nurdin H I and Yamamoto N 2020 Phys. Rev. Applied 14(2) 024065
  * [30] Pfeffer P, Heyder F and Schumacher J 2022 Phys. Rev. Research 4(3) 033176
  * [31] Mujal P, Martínez-Peña R, Giorgi G L, Soriano M C and Zambrini R 2023 npj Quantum Inf. 9 16 ISSN 2056-6387 URL https://doi.org/10.1038/s41534-023-00682-z
  * [32] Čindrak S, Donvil B, Lüdge K and Jaurigue L C 2024 Phys. Rev. Res. 6(1) 013051 URL https://link.aps.org/doi/10.1103/PhysRevResearch.6.013051
  * [33] Hülser T, Köster F, Jaurigue L C and Lüdge K 2022 Opt. Mater. Express 12 1214
  * [34] Huang G B, Zhu Q Y and Siew C K 2004 IEEE proceedings 2 985–990
  * [35] Biasi S, Franchi R, Cerini L and Pavesi L 2023 APL Photonics 8 096105 ISSN 2378-0967 (Preprint "") URL https://doi.org/10.1063/5.0156189
  * [36] Lorenz E N 1963 J. Atmos. Sci. 20 130
  * [37] Kantz H and Schreiber T 2003 Nonlinear Time Series Analysis 2nd ed (Cambridge University Press)
  * [38] Tsuchiyama K, Röhm A, Mihana T, Horisaki R and Naruse M 2023 Chaos 33 063145 ISSN 1054-1500 URL https://doi.org/10.1063/5.0143846
  * [39] Storm L, Gustavsson K and Mehlig B 2022 Mach. Learn.: Sci. Technol. 3 045021
  * [40] Viehweg J, Worthmann K and Mäder P 2023 Neurocomputing 522 214–228 ISSN 0925-2312 URL https://www.sciencedirect.com/science/article/pii/S0925231222014291
  * [41] Lukosevicius M and Jaeger H 2009 Computer Science Review 3 127–149 ISSN 1574-0137
  * [42] Hart A G, Hook J L and Dawes J H 2020 Neural Netw. 128 234–247 ISSN 0893-6080
  * [43] Jaurigue L and Lüdge K 2024 Neuromorph. Comput. Eng. 4 014001 URL https://dx.doi.org/10.1088/2634-4386/ad1d32


Recurrent neural networks are capable of learning the dynamics of an unknown nonlinear system purely from input-output measurements. However, the resulting models do not provide any stability guarantees on the input-output mapping. In this work, we represent a recurrent neural network as a linear time-invariant system with nonlinear disturbances. By introducing constraints on the parameters, we can guarantee finite gain stability and incremental finite gain stability. We apply this identification method to learn the motion of a four-degrees-of-freedom ship that is moving in open water and compare it against other purely learning-based approaches with unconstrained parameters. Our analysis shows that the constrained recurrent neural network has a lower prediction accuracy on the test set, but it achieves comparable results on an out-of-distribution set and respects stability conditions.

Recurrent Neural Networks, System Identification, Linear Matrix Inequality Constraints, Deep Learning

§ INTRODUCTION

§ RELATED WORK

§ BACKGROUND

§ METHODOLOGY

§ NUMERICAL EXPERIMENT

§ CONCLUSION

We thank Alexandra Baier for helpful discussions and her python package deepsysid that was used for training, testing, and evaluating our models. This work is funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2075 – 390740016. We acknowledge the support by the Stuttgart Center for Simulation Science (SimTech). The authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS)
for supporting Daniel Frank.

[Baier and Staab, 2022]
Alexandra Baier and Steffen Staab.
A Simulated 4-DOF Ship Motion Dataset for System Identification
under Environmental Disturbances, 2022.
URL <https://doi.org/10.18419/darus-2905>.

[Baier et al., 2021]
Alexandra Baier, Zeyd Boukhers, and Steffen Staab.
Hybrid Physics and Deep Learning Model for Interpretable Vehicle
State Prediction.
arXiv preprint arXiv:2103.06727, 2021.

[Boyd and Vandenberghe, 2004]
Stephen Boyd and Lieven Vandenberghe.
Convex Optimization.
Cambridge University Press, 2004.

[Fazlyab et al., 2019]
Mahyar Fazlyab, Alexander Robey, Hamed Hassani, Manfred Morari, and George J.
Efficient and Accurate Estimation of Lipschitz Constants for Deep
Neural Networks.
33rd Conference on Neural Information Processing Systems, 32,

[Goodfellow et al., 2016]
Ian Goodfellow, Yoshua Bengio, and Aaron Courville.
Deep Learning.
MIT Press, 2016.

[Gu et al., 2022]
Fangda Gu, He Yin, Laurent El Ghaoui, Murat Arcak, Peter Seiler, and Ming Jin.
Recurrent Neural Network Controllers Synthesis with Stability
Guarantees for Partially Observed Systems.
Proceedings of the AAAI Conference on Artificial Intelligence,
360 (5):0 5385–5394, 2022.

[Hasselmann et al., 1973]
K. Hasselmann, T. P. Barnett, E. Bouws, H. Carlson, D. E. Cartwright, K. Enke,
J. A. Ewing, H. Gienapp, D. E. Hasselmann, P. Kruseman, A Meerburg,
A. Müller, D. J. Olbers, K. Richter, W. Sell, and H. Walden.
Measurements of Wind-Wave Growth and Swell Decay during the Joint
North Sea Wave Project (JONSWAP).
Ergänzungsheft zur Deutschen Hydrographischen Zeitschrift,
Reihe A, 1973.

[Hochreiter and Schmidhuber, 1997]
Sepp Hochreiter and Jürgen Schmidhuber.
Long Short-Term Memory.
Neural computation, 90 (8):0 1735–1780, 1997.

[Isherwood, 1972]
R.M. Isherwood.
Wind Resistance of Merchant Ships.
The Royal Institution of Naval Architects, 115:0
327–338, 1972.

[Junnarkar et al., 2022]
Neelay Junnarkar, He Yin, Fangda Gu, Murat Arcak, and Peter Seiler.
Synthesis of Stabilizing Recurrent Equilibrium Network Controllers.
arXiv preprint arXiv:2204.00122, 2022.

[Kingma and Ba, 2015]
Diederik P. Kingma and Jimmy Lei Ba.
Adam: A Method for Stochastic Optimization.
In 3rd International Conference on Learning Representations
(Poster), 2015.

[Mohajerin and Waslander, 2019]
Nima Mohajerin and Steven L. Waslander.
Multistep Prediction of Dynamic Systems With Recurrent Neural
IEEE Transactions on Neural Networks and Learning Systems,
300 (11):0 3370–3383, 2019.

[Pauli et al., 2022]
Patricia Pauli, Julian Berberich, and Frank Allgöwer.
Robustness analysis and training of recurrent neural networks using
dissipativity theory.
at-Automatisierungstechnik, 700 (8):0
730–739, 2022.

[Perez et al., 2006]
Tristan Perez, Andrew Ross, and Thor Fossen.
A 4-DOF simulink model of a coastal patrol vessel for manoeuvring in
In Proceedings of the 7th IFAC Conference on Manoeuvring and
Control of Marine Craft, pages 1–6. International Federation for Automatic
Control, 2006.

[Pillonetto et al., 2022]
Gianluigi Pillonetto, Tianshi Chen, Alessandro Chiuso, Giuseppe De Nicolao, and
Lennart Ljung.
Classical System Identification.
In Regularized System Identification, pages 17–31. Springer,

[Revay et al., 2021]
Max Revay, Ruigang Wang, and Ian R. Manchester.
A Convex Parameterization of Robust Recurrent Neural Networks.
IEEE Control Systems Letters, 50 (4):0
1363–1368, 2021a.

[Revay et al., 2021]
Max Revay, Ruigang Wang, and Ian R. Manchester.
Recurrent equilibrium networks: Unconstrained learning of stable and
robust dynamical models.
In 2021 60th IEEE Conference on Decision and Control (CDC),
Austin, TX, USA, December 14-17, 2021, pages 2282–2287. IEEE,
URL <https://doi.org/10.1109/CDC45484.2021.9683054>.

[Sastry, 2013]
Shankar Sastry.
Nonlinear systems: analysis, stability, and control,
volume 10.
Springer Science & Business Media, 2013.

[Scherer and Weiland, 2000]
Carsten Scherer and Siep Weiland.
Linear matrix inequalities in control.
Lecture Notes, Dutch Institute for Systems and Control, Delft,
The Netherlands, 30 (2), 2000.

[Scherer, 2006]
Carsten W. Scherer.
LMI Relaxations in Robust Control.
European Journal of Control, 120 (1):0 3–29,

[Suykens et al., 1995]
Johan A. K. Suykens, Bart L. R. De Moor, and Joos Vandewalle.
Nonlinear system identification using neural state space models,
applicable to robust control design.
International Journal of Control, 620 (1):0
129–152, 1995.

[Veenman et al., 2016]
Joost Veenman, Carsten W. Scherer, and Hakan Köroglu.
Robust stability and performance analysis based on integral quadratic
European Journal of Control, 31:0 1–32, 2016.
URL <https://doi.org/10.1016/j.ejcon.2016.04.004>.

§ HYPERPARAMETER OPTIMIZATION