# XRJL-HKUST at SemEval-2021 Task 4: WordNet-Enhanced Dual Multi-head Co-
Attention for Reading Comprehension of Abstract Meaning

Yuxin Jiang*, Ziyi Shou*, Qijun Wang, Hao Wu, Fangzhen Lin
HKUST-Xiaoi Robot Joint Lab on Machine Learning and Cognitive Reasoning
Department of Computer Science and Engineering
The Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong
{yjiangcm,zshou,qwangcd,hwubx}<EMAIL_ADDRESS>

###### Abstract

This paper presents our submitted system to SemEval 2021 Task 4: Reading
Comprehension of Abstract Meaning. Our system uses a large pre-trained
language model as the encoder and an additional dual multi-head co-attention
layer to strengthen the relationship between passages and question-answer
pairs, following the current state-of-the-art model DUMA. The main difference
is that we stack the passage-question and question-passage attention modules
instead of calculating parallelly to simulate re-considering process. We also
add a layer normalization module to improve the performance of our model.
Furthermore, to incorporate our known knowledge about abstract concepts, we
retrieve the definitions of candidate answers from WordNet and feed them to
the model as extra inputs. Our system, called WordNet-enhanced DUal Multi-head
Co-Attention (WN-DUMA), achieves 86.67% and 89.99% accuracy on the official
blind test set of subtask 1 and subtask 2 respectively.

##  1 Introduction

Recently, there has been an increasing interest on Machine Reading
Comprehension (MRC).††*Equal contribution. While most MRC studies such as
CNN/Daily Mail Hermann et al. (2015) focus on concrete concepts, SemEval 2021
Task 4, Reading Comprehension on Abstract Meaning (ReCAM), targets abstract
concept understanding, including imperceptibility in subtask 1 and
nonspecificity in subtask 2. The former, imperceptibility, highlights the
abstract words that refer to ideas and concepts that do not correspond
directly to human perception. The latter is for hypernyms and abstract
concepts such as the class of vertebrate which includes whales as a concrete
subclass Changizi (2008).

In this task, given news fragments and incomplete abstracts, the machine needs
to select the most suitable abstract words from candidate answers. Figure 1
shows one example of ReCAM subtask 1.

Figure 1: An example of ReCAM subtask 1.

Passage is the news sections and Question is a human written summary in which
abstract words have been removed. Machines are requested to choose abstract
words from five candidates for replacing @placeholder.

For this shared task, we regard both subtasks as multi-choice MRC tasks.
Various deep neural networks and attention mechanisms (e.g. Dhingra et al.
(2017); Wang et al. (2018); Zhang et al. (2020a, b); Jin et al. (2020)) have
been proposed to address these tasks. In our work, following the state-of-the-
art model DUMA Zhu et al. (2020), we adopt a Pre-trained Language Model (PrLM)
as encoder and extend with an additional dual multi-head co-attention layer to
strengthen the relationship between passages and question-answer pairs. For
the dual multi-head attention layer, while DUMA builds passage-question and
question-passage attention modules in a parallel way to simulate the
transposition thinking process, our model stacks two attention modules in
order to simulate the process of re-considering for a deeper understanding of
the passage. More details on our attention calculation process can be found in
Section 2. Furthermore, we add an additional layer normalization module
immediately after the attention module. From our experiments, we found that
this additional normalization module definitely improves our model’s
performance.

Our most significant design decision is to use WordNet Miller (1995) to expand
on the abstract concepts in the candidate answers. Intuitively, expanding an
abstract concept according to its definition in a dictionary should help as it
helps relate the abstract concept with others that may occur in the text. A
key conclusion that we can draw from our experiments is that this is indeed
the case. One problem that we encountered when implementing this idea was that
most English words have multiple entries in WordNet. For example, bank in
WordNet can have Noun definitions as well as Verb definitions. We addressed
this problem by using some heuristics and some additional information such as
part-of-speech labels. Because of the significant role played by WordNet, we
call our system WordNet-enhanced DUal Multi-head Co-Attention (WN-DUMA).

We remark that our system did not use any additional training data for the
tasks. In the final evaluation, our model is ranked 10 out of 23 and 9 out of
28 on the official subtask 1 and subtask 2 blind test set with 86.67% and
89.99% accuracy, respectively. The code for our model is publicly
available111https://github.com/zzshou/RCAM.

The rest of the paper is organized as follows. Section 2 gives the details of
our system. Section 3 describes our experimental setup including the datasets
and hyper parameters used for training. Section 4 presents experimental
results. Section 5 concludes this paper with some final remarks.

##  2 System Description

In this section, We describe the framework of our end-to-end model WN-DUMA.
Figure LABEL:model_architecture depicts the detailed architecture of our
approach, with inputs at the bottom and outputs at the top.

#### WordNet-enhanced Encoder

We regard both subtask 1 and 2 as multi-choice MRC problems. Such a problem
includes a passage, a question with a @placeholder, and 5 candidate answers to
choose from. First, we replace @placeholder in the question with the given 5
candidate answers to form 5 options. In the tasks, the candidate answers are
all single words with abstract meanings, so we decided to add some extra
knowledge from WordNet Miller (1995) to help the system better understanding
the abstract meanings. More specifically, for a single candidate answer, we
find its part-of-speech tag based on the option it’s located in, and extract
its definitions under this part-of-speech tag. After tokenization, every
instance is cast into the input form: [CLS] _passage_ [SEP] _option + answer
definition_ [SEP]. To encode input tokens into representations, we feed them
through a PrLM based on Transformer to obtain sequence embeddings, which draws
a global relationship between the passage and the option-definition.

#### Dual Multi-head Co-Attention Layer

Based on the above process, we further separate the output representations
from transformer encoder to acquire the passage context embeddings
$E^{P}\in\mathbb{R}^{d_{model}\times l_{p}}$ and the context embeddings of
option-definition $E^{OD}\in\mathbb{R}^{d_{model}\times l_{od}}$, where
$l_{p},l_{od}$ denote the maximum length of passage and option-definition
respectively. Then based on the bi-directional matching network of DUMA which
is quite similar to the multi-head self-attention module in vanilla
transformer block Vaswani et al. (2017), we first take $E^{OD}$ as Query,
$E^{P}$ as Key and Value to calculate one of the co-attention representations,
which simulates the process of human re-reading the passage with impression of
option and definition. The formulas are listed as follows:

$\displaystyle Q_{i}=E^{OD}W_{i}^{Q}$ (1) $\displaystyle K_{i}=E^{P}W_{i}^{K}$
(2) $\displaystyle V_{i}=E^{P}W_{i}^{V}$ (3) $\displaystyle
head_{i}=softmax(\frac{Q_{i}K_{i}^{T}}{\sqrt{d_{k}}})V_{i}$ (4) $\displaystyle
MHA=Concat(head_{1},...,head_{h})W^{O}$ (5) $\displaystyle
REP_{1}=Normalize(E^{OD}+MHA)$ (6)

where $W_{i}^{Q}\in\mathbb{R}^{d_{model}\times d_{q}}$,
$W_{i}^{K}\in\mathbb{R}^{d_{model}\times d_{k}}$,
$W_{i}^{V}\in\mathbb{R}^{d_{model}\times d_{v}}$,
$W_{i}^{O}\in\mathbb{R}^{hd_{v}\times d_{model}}$ are linear transformations
with learnable parameters, $d_{q},d_{k},d_{v}$ denote the dimension of _Query_
, _Key_ and _Value_ , $h$ denotes the number of heads. Different from the
structure of DUMA, here we make two changes: 1) apply ”Add and Normalize”
after getting the multi-head attention representation, which could result in
more stable training. 2) compute another co-attention representation by
stacking rather than paralleling: take the acquired $REP_{1}$ as _Key_ and
_Value_ , $E^{P}$ as _Query_ , which simulates the process of re-considering
the option-definition with deeper understanding of the passage. Finally, we
obtain $REP_{1}$ and $REP_{2}$, which have the same dimension as $E^{OD}$ and
$E^{P}$, respectively. As a result, we can stack the co-attention module for k
layers.

#### Classifier

Here the co-attention representations $REP_{1}$ and $REP_{2}$ are merged and
used for final classification:

$\displaystyle I_{1}=AvgPool(REP_{1})$ (7) $\displaystyle
I_{2}=AvgPool(REP_{2})$ (8) $\displaystyle M=Concat(I_{1},I_{2})$ (9)
$\displaystyle logits=MW_{M}$ (10)

where $I_{1},I_{2}\in\mathbb{R}^{d_{model}}$, $M\in\mathbb{R}^{2d_{model}}$,
$W_{M}\in\mathbb{R}^{d_{model}\times n_{class}}$ denotes the one-layer fully-
connected neural network, $n_{class}$ denotes the number of candidate answers.
Consequently, for a single instance, we could get as many logits as the
candidate answers, which are used to compute the corss-entropy loss by
softmax.

##  3 Experimental Setup

#### Data and Metric

We used the official datasets Zheng et al. (2021) provided by SemEval 2021
Task 4 competition. They were collected from BBC News in English language.
Some basic statistics are listed in Table 1.

| Task 1 | Task 2
---|---|---
Train | 3,227 | 3,318
Dev | 837 | 851
Test | 2,025 | 2,017
Avg. passage length | 270.3 | 429.7
Avg. question length | 24.6 | 27.1
Vocabulary size | 16,318 | 17,006
Answer vocabulary size | 4,333 | 4,775
Table 1: Basic statistics of subtask 1 and subtask 2 dataset.

According to the requirement of the organizers, participants could only use
the corresponding dataset for a specific subtask to build models to ensure
fairness. For better performance, technics like multi-task learning Wan (2020)
are recommended for MRC tasks. In both subtask 1 and subtask 2, we utilize
accuracy as the metric to evaluate our model performance.

#### Hyper Parameters

All of our codes are written based on PyTorch222https://pytorch.org/. To
extract the word definition of candidate answers, we use NLTK toolkit Bird et
al. (2009). The transformer encoder we used is pretrained ALBERT-xxlarge-v2
model333https://github.com/huggingface/transformers. Since the code of DUMA is
not open-source, we reimplement it by only using one co-attention layer where
the attention heads are 64 and the dimension of _Query_ , _Key_ and _Value_
are all 64, because it is pointed that more co-attention layers do not improve
the performance Zhu et al. (2020). The setting is also applied to our WN-DUMA
for fair comparison.

Due to limited resources, the maximum sequence length of input tokens is set
to 150 for both subtask 1 and subtask 2. In fact, we found that sequence
length longer than 150 can only slightly improve the model performance. We
choose mini-batch size equal to 2, and the AdamW optimizer Loshchilov and
Hutter (2018) with an initial learning rate of 5e-06. We use some strategies
for more stable training: 1) clip the gradient norm to 10; 2) adopt a linear
scheduler with warm up of the first 10% training steps. To avoid overfitting,
we apply 0.1 dropout Srivastava et al. (2014) rate to the co-attention layer.
We trained all the models for 3 epochs, evaluate on the dev set at every 200
training steps and save the model with the best dev accuracy. For each single
model, we run experiments for 5 times with different random seeds and use the
average as the ultimate performance.

##  4 Results

Model | Task 1 dev | Task 1 test | Task 2 dev | Task 2 test
---|---|---|---|---
$\rm BERT_{large}$ Devlin et al. (2019) | 67.74 | - | 69.45 | -
$\rm RoBERTa_{large}$ Liu et al. (2019) | 74.31 | - | 74.50 | -
$\rm ALBERT_{xxlarge}$ Lan et al. (2020) | 84.83 | - | 82.84 | -
$\rm ALBERT_{xxlarge}$ \+ DUMA Zhu et al. (2020) | 85.07 | - | 86.13 | -
$\rm ALBERT_{xxlarge}$ (question only) | 79.57 | - | 82.14 | -
$\rm ALBERT_{xxlarge}$ \+ WN-DUMA (single) | 85.90 | 84.54 | 87.43 | 86.61
$\rm ALBERT_{xxlarge}$ \+ WN-DUMA (ensemble) | 87.57 | 86.67 | 90.01 | 89.99
Table 2: Model comparison on subtask 1 and subtask 2 dataset.

###  4.1 Quantitative Analysis

Table 2 summarizes the experimental results. The first three models only have
encoder (without enhancement of WordNet) and classifier part. It is clearly
seen that ALBERT is much more efficient as encoder for abstract meaning
understanding. It is worth noting that by only using the question and answer
as input, the ALBERT model can also get pretty good results, as table 2 shows.
Intuitively, it may be because the model could utilize syntax and semantics of
the question sentence to choose the correct answer without looking through the
passage.

Compared to $\rm ALBERT_{xxlarge}$, adding DUMA layer obtains around 0.2%
improvement in subtask 1, and more than 3% improvement in subtask 2. Besides,
our WN-DUMA single model achieves further improvements based on DUMA on both
subtasks, +0.83% and +1.3% respectively, without increasing the number of
parameters. Using a majority vote scheme, we ensemble our WN-DUMA model with
different parameters for more stable predictions. Eventually, our ensemble
models which get 87.57% on subtask 1 dev set and 90.01% on subtask 2 dev set
acquire the best performance on test sets (86.67% and 89.99%, respectively)
among our submissions.

Figure 3 and Figure 4 illustrate the dev accuracy of different models on
subtask 1 and subtask 2 as the number of training steps increases. It is
interesting to observe that models with co-attention layer (DUMA and WN-DUMA)
could get over 70% accuracy with only 10% of training examples. While ALBERT
model has to be trained with the full dataset to get relatively high accuracy.
Consequently, our WN-DUMA model may be useful when there only exists a small
amount of training data.

Figure 3: Subtask 1 dev accuracy over number of training steps. Figure 4:
Subtask 2 dev accuracy over number of training steps.

###  4.2 Error Analysis

In order to further improve our model performance in the future, we analyze
some incorrect predictions made by WN-DUMA, and classify them into two
categories:

  * •

Candidate answers with similar meanings. In some failure cases, the
similarities between candidate answers are too high to distinguish. For
example, outstanding and extraordinary, challenge and attempt, etc.

  * •

Lack of commonsense and relying too much on the information of the passage.
Due to the fact that the question is the summary of the passage, the machine
need to choose the most appropriate answer from a global perspective with some
commonsense. However, our model make decisions by only capturing the local
information in some cases. A specific example can be seen in Figure 5. We can
see that the model predicts that the answer is ”troubled”, most likely because
the passage mentions ”the school was trapped into financial difficulties”.

Figure 5: An failure example made by WN-DUMA. The ground true is the answer
with a correct mark at the bottom. While the prediction is the answer with
”predicted” at its bottom.

##  5 Conclusion

In this paper, we describe our submitted system in SemEval 2021 Task 4 ReCAM.
Unlike previous MRC datasets, ReCAM focus more on machine’s ability in
understanding and representing abstract concepts. In order to provide more
knowledge of abstract word, we extract WordNet definitions for each candidate
answer based on part-of-speech tags. In addition, our proposed WN-DUMA model
consists of a PrLM as the encoder and a dual multi-head co-attention layer to
enhance the relationship between passage and question-answer pairs as human’s
re-considering process. Our WN-DUMA model improves the performance of our
baseline model DUMA on these datasets.

There are some limitations in our experiments. Firstly, training data size of
this task is limited compared to other MRC tasks, with less than 3400 training
pairs in both subtasks. This is understandable as collecting labeled data in
many natural language processing tasks is expensive. Secondly, using $\rm
ALBERT_{xxlarge}$ PrLM, we only set 150 as the maximum text length in our
experiments due to device limitation. Important sentences in the passage that
are highly relevant to the summary are sometimes not covered. For PrLMs, their
performance always improve as the number of their parameters increase. The use
of large pre-trained models sometimes requires the sacrifice of context. For
our future work, we plan to explore ways to train models more efficiently with
limited amount of labeled data, and to design more cost-effective models to
deal with long input texts.

## References

  * Bird et al. (2009) Steven Bird, Ewan Klein, and Edward Loper. 2009.  _Natural Language Processing with Python_.  O’Reilly Media.
  * Changizi (2008) Mark A Changizi. 2008.  Economically organized hierarchies in wordnet and the oxford english dictionary.  _Cognitive Systems Research_ , 9(3):214–228.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Dhingra et al. (2017) Bhuwan Dhingra, Hanxiao Liu, Zhilin Yang, William W. Cohen, and Ruslan Salakhutdinov. 2017.  Gated-attention readers for text comprehension.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics, ACL 2017, Vancouver, Canada, July 30 - August 4, Volume 1: Long Papers_ , pages 1832–1846. Association for Computational Linguistics.
  * Hermann et al. (2015) Karl Moritz Hermann, Tomáš Kočiskỳ, Edward Grefenstette, Lasse Espeholt, Will Kay, Mustafa Suleyman, and Phil Blunsom. 2015.  Teaching machines to read and comprehend.  In _Proceedings of the 28th International Conference on Neural Information Processing Systems-Volume 1_ , pages 1693–1701.
  * Jin et al. (2020) Di Jin, Shuyang Gao, Jiun-Yu Kao, Tagyoung Chung, and Dilek Hakkani-tur. 2020.  Mmm: Multi-stage multi-task learning for multi-choice reading comprehension.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34(05):8010–8017.
  * Lan et al. (2020) Zhenzhong Lan, Mingda Chen, Sebastian Goodman, Kevin Gimpel, Piyush Sharma, and Radu Soricut. 2020.  ALBERT: A lite BERT for self-supervised learning of language representations.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_.
  * Loshchilov and Hutter (2018) Ilya Loshchilov and Frank Hutter. 2018.  Fixing weight decay regularization in adam.
  * Miller (1995) George A. Miller. 1995.  Wordnet: A lexical database for english.  _Commun. ACM_ , 38(11):39–41.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey E. Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. 2014.  Dropout: a simple way to prevent neural networks from overfitting.  _J. Mach. Learn. Res._ , 15(1):1929–1958.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , pages 5998–6008.
  * Wan (2020) Hui Wan. 2020.  Multi-task learning with multi-head attention for multi-choice reading comprehension.  _arXiv preprint arXiv:2003.04992_.
  * Wang et al. (2018) Shuohang Wang, Mo Yu, Jing Jiang, and Shiyu Chang. 2018.  A co-matching model for multi-choice reading comprehension.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics, ACL 2018, Melbourne, Australia, July 15-20, 2018, Volume 2: Short Papers_ , pages 746–751. Association for Computational Linguistics.
  * Zhang et al. (2020a) Shuailiang Zhang, Hai Zhao, Yuwei Wu, Zhuosheng Zhang, Xi Zhou, and Xiang Zhou. 2020a.  Dcmn+: Dual co-matching network for multi-choice reading comprehension.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34(05):9563–9570.
  * Zhang et al. (2020b) Zhuosheng Zhang, Yuwei Wu, Junru Zhou, Sufeng Duan, Hai Zhao, and Rui Wang. 2020b.  Sg-net: Syntax-guided machine reading comprehension.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34(05):9636–9643.
  * Zheng et al. (2021) Boyuan Zheng, Xiaoyu Yang, Yuping Ruan, Quan Liu, Zhen-Hua Ling, Si Wei, and Xiaodan Zhu. 2021.  SemEval-2021 task 4: Reading comprehension of abstract meaning.  In _Proceedings of the 15th International Workshop on Semantic Evaluation (SemEval-2021)_.
  * Zhu et al. (2020) Pengfei Zhu, Hai Zhao, and Xiaoguang Li. 2020.  Duma: Reading comprehension with transposition thinking.  _arXiv preprint arXiv:2001.09415_.


# Backpropagation on Dynamical Networks

Eugene Tan, Dèbora Corrêa, Thomas Stemler, Michael Small E. Tan, T. Stemler
and M. Small are with the Department of Mathematics & Statistics, The
University of Western Australia, Crawley, WA 6009, Australia
E-mail<EMAIL_ADDRESS>D. Corrêa is with the Department of Computer
Science & Software Engineering, The University of Western Australia, Crawley,
WA 6009.Manuscript received INSERT DATE; revised INSERT DATE.

###### Abstract

Dynamical networks are versatile models that describe a variety of behaviours
such as synchronisation and feedback. However, applying these models in real
systems is difficult as prior information of the connectivity structure or
local dynamics is often unknown and must be inferred from observations of
network states. Additionally, the influence of coupling interactions further
complicates the isolation of local node dynamics. Given the architectural
similarities between dynamical networks and recurrent neural network (RNNs),
we propose a network inference method based on the backpropagation through
time algorithm used to train RNNs. This method aims to simultaneously infer
both the connectivity structure and local node dynamics from node state
observations. An approximation of local node dynamics is first constructed
using a neural network. This is alternated with an adapted BPTT algorithm to
regress corresponding network weights by minimising prediction errors of the
network based on the previously constructed local models until convergence.
This method was found to be successful in identifying the connectivity
structure for coupled networks of chaotic oscillators. Freerun prediction
performance with the resulting local models and weights was comparable to the
true system with noisy initial conditions. The method is also extended to
asymmetric negative coupling.

###### Index Terms:

dynamical networks, network inference, backpropagation, neural networks,
machine learning

##  1 Introduction

Dynamical networks are a common occurrence when studying complex systems where
the goal is to model large systems of multiple interacting components. In its
simplest form, a dynamical network may be defined as having three main
components: local node dynamics, coupling dynamics and connectivity structure
(see Equation 1). One key feature of dynamical networks is its ability to
recreate and describe rich and interesting dynamics such as chimera states,
synchronisation [1, 2, 3] and cascading failure commonly encountered in real
world systems [4]. They have also been applied to other systems such as neuron
networks [5, 4], power grids [6, 7], epidemic spread [8], and cardiac
arrhythmia [9].

$\dot{x}_{i}=f(x_{i})+\sum_{i\neq j}c_{i,j}g(x_{i},x_{j}).$ (1)

In Eq 1, nodes are indexed with $i$, $c$ represents the coupling structure,
$f$ and $g$ represents the local and coupling dynamics respectively. Despite
its relatively simple form, data-driven applications of this framework to
model real world systems can be difficult if information regarding the
coupling structure and dynamics is not accessible. In these cases, these
components need to be inferred before a dynamical networks framework can be
applied. One simple example would be the data driven analysis of a network of
chaotic oscillators. Specifically, there is the inverse problem of inferring
network characteristics from data. However, this problem is difficult due to
the complicated interplay between local models and coupling behaviour. Many
existing methods of inferring dynamical network properties require either the
local or coupling dynamics to be at least partially known. Alternatively, some
methods restrict inferences to statistical arguments with uncertainty in order
to remove the dependence on a priori information.

Numerous approximation and inference methods for inferring network
connectivity have been developed each with varying success. These methods may
be classified into three main approaches:

  1. 1.

Direct Approach \- For systems where both local $f$ and coupling dynamics $g$
is explicitly known, a direct algebraic approach using derivatives have been
used to recover connectivity structure [10].

  2. 2.

Perturbation Methods \- These methods try and recreate a proxy dynamical
network that performs similarly to the target system. Small perturbations are
applied to individual nodes in the network and observed to see how each signal
propagates. These propagations are then used to infer the network connectivity
structure [11, 12, 13] .

  3. 3.

Signal Correlation \- This approach focuses on quantifying the correlation
between node signals [14]. Node pairs with highly correlated signals are used
to infer the connectivity structure. Maximum route entropy has also been used
as an alternative to achieve the same result. [15]. An extension of this
approach is the employment of some notion of causality such as Granger
causality [16, 17] and short term causal dependence. [11].

The inference of connectivity structure also poses an additional challenge of
false positive couplings and remains an open research area. Granger causality
and similar statistical based approaches have been found to particularly
struggle with this problem, even for simple network structures such as rings
and chains [17, 18].

The second problem of constructing models to describe the local dynamics $f$
in multi-node networks has been studied for pseudo-periodic systems [9, 16].
Outside these, many approaches focus on the case of a single node without
external interactions and utilise a time-delay embedding with appropriately
chosen delay and embedding dimension [19]. Proposed modelling methods for long
term prediction include neural network predictors [20], radial basis
predictors [21] and reservoir computing [22] among several others. Assessing
the quality of these local models relies on the comparison of dynamical
invariants such as Lyapunov exponents or correlation dimension [23, 24, 25]
and attractor homology [26].

There is also a third related problem of identifying coupling equations $g$
from signals. Solutions have been proposed using phase dynamics for networks
exhibiting pseudo-periodic behaviour, but not necessarily for aperiodic
behaviour [9]. For simple oscillator networks, Panaggio et al. [27] propose a
gradient descent method to regress coupling functions of a given form. This
approach is draws similarities with [10]. For simplicity, this paper assumes
that the coupling equations are known.

For low dimensional systems, symbolic approaches aimed at tackling either the
dynamics or connectivity structure inference problem has gained increased
traction in recent times. A method based on the ideas of compressed sensing
[28] and sparse optimisation methods (SINDy) [29] have been proposed to
reconstruct the functional form of the dynamics equations using a library of
basis functions (e.g. power series, Fourier series) for isolated dynamical
systems with no interactions. Compressive sensing tries to find the sparsest
combination of coefficients for some power series basis via convex
optimisation with respect to the $L_{1}$-norm. Link and coupling weight
identification can also be included into this approach [28, 30]. However, we
note that this method may be limited to simple systems with analytical or well
behaved equations.

The symbolic approach with basis functions was also extended to local dynamics
and coupling functions by utilising prior known information of the network
connectivity structure [31]. A similar approach (ARNI) [32] using basis
functions has also been proposed to reconstruct the coupling function and
connectivity structure of dynamical networks. ARNI is also applicable to
hetereogeneous networks where local dynamics $f$ differ between nodes. We note
that ARNI can also be adapted to infer local dynamics instead [31].

Despite these topics being well studied, the literature on inferring both
connectivity structure and node dynamics from time series observations is not
rich. Methods for reconstructing local dynamics (e.g. SINDy [29] and Gao [31])
are either restricted to single node systems, or rely or prior information
about the connectivity structure. In cases where the aim is to infer
connectivity structure, either local dynamics are known prior [10] or are not
directly inferred [32]. There remains the need of a method that tackles the
inference of both of these properties simultaneously.

We note that a method utilising a mean field approach has been proposed to
tackle the simultaneous inference of local dynamics and connectivity structure
[33]. This method constructs an ‘effective network’ that acts as an
approximate proxy for the real dynamical network, reproducing many of its
properties and behaviours. Under specific assumptions such as weak coupling,
network heterogeneity and limited links, this method has been shown to perform
well in identifying critical transitions and the connectivity structure of a
cat cortex. However, the ‘effective network’ relies on the network sparsity
and presence of hubs to recover an approximation of the local dynamics.
Additionally, this approximation of local dynamics is not fully disentangled
from the interaction dynamics.

In this paper, we present a proof of principle of a data-driven approach for
simultaneously inferring both the local node dynamics and weighted
connectivity structure from the observation of node signals, assuming only
prior knowledge of the coupling function. Recognising the structural
similarities between dynamical networks and RNNs, we propose an adaptation of
the backpropagation through time (BPTT) algorithm used to train RNNs [34], as
a method of regressing the connectivity structure of a dynamical network. This
regression is then used to decouple the observed node signals and
simultaneously construct a model of local dynamics. This method combines the
tasks achieved by several of the abovementioned methods, whilst also dealing
instances of dynamical networks with weighted edges.

We test the performance of the proposed method on a dynamical network of
Lorenz [35] and Chua oscillators [36] with a randomly initialised weighted
connectivity structure. We find that an adapted backpropagation regression
method was able to reproduce the network coupling structure and local dynamics
for a 64 node network. Additional cases with asymmetric and negatively
weighted coupling are also explored, yielding similar results. To test the
performance on a more realistic application, the proposed backpropagation
regression method was applied to a network of driven FitzHugh-Nagumo
oscillators operating in the chaotic regime to simulate a simple biological
neuron network. Similar to the previous cases, the backpropagation method was
able to successfully reproduce both the network structure and local dynamics.

##  2 Background

###  2.1 Dynamical Networks

Dynamical networks can be defined as a network $G=(f,g,C)$ consisting of three
main components, local dynamics $f$, coupling dynamics $g$ and connectivity
structure $C$. Node states evolve according to Equation 1 with local behaviour
governed by $f$ and additional external influences from neighbours due to
coupling components $g$ and $C$.

The effect of coupling between nodes results in a high dimensional system with
multiple interacting sub-systems. For small or highly connected networks,
sufficiently large coupling can cause network-wide synchronisation of node
states. As coupling is decreased, chimera states and group synchronisation
emerges [1, 2, 3]. When no coupling is present, all nodes behave independent
of each other.

Here, we focus on dynamical networks based on two different 3-dimensional
chaotic systems, Lorenz and Chua, with diffusive coupling in the first
component. For the latter, we will focus on the Chua system with cubic
nonlinearity $\phi$ operating in the single scroll mode [36] to analyse
dynamics contrasting with the Lorenz system.

###  2.2 Backpropagation and Recurrent Neural Networks

Recurrent neural networks (RNNs) are a modelling architecture developed in the
field of machine learning that is widely used for sequence and time series
prediction [37, 38, 39]. RNNs consist of an input layer with weights
$C^{(in)}$, a fully connected hidden layer of nodes $\mathbf{h}(t)$ with
feedback weights $C^{(h)}$, and an output layer with readout weights
$C^{(out)}$. For prediction, an input sequence $\mathbf{x}(t_{n})$ is fed into
the RNN via the input layer with input weights $C^{(in)}$. These inputs are
usually discretely sampled at regular time intervals ${t_{n}}$ and then
sequentially fed into the recurrent hidden layer whose nodes feedback to each
other according to weights before applying a nonlinear transformation
$\sigma$. In its continuous time formulation, this allows network to
approximate nonlinear dependencies between states with nodes evolving
according to Equation 2,

$\mathbf{\hat{h}}(t)=\sigma(C^{(in)}\mathbf{x}(t)+C^{(h)}\mathbf{h}(t)+\mathbf{b}).$
(2)

Predictions $\hat{\mathbf{x}}(t)$ are calculated using the readout weights and
hidden node states as given in Equation 3,

$\mathbf{\hat{x}}(t)=C^{(out)}\cdot\mathbf{h}(t).$ (3)

The training of RNNs involve the small adjustment of network weights
$C^{(in)}$, $C^{(h)}$, $C^{(out)}$ and biases $\mathbf{b}$. This is achieved
via the backpropagation through time (BPTT) algorithm [16]. The BPTT algorithm
can be interpreted as an adaptation of the gradient descent training algorithm
for neural networks that accounts for the propagation of network weight
effects through the history of the RNN’s hidden node states.

The BPTT algorithm works by unfolding the predictions of the recurrent network
at each time step. This process reframes the forward pass of an RNN into a
feedforward network where each layer corresponds to the node states of the RNN
at a specific time. For a given RNN model, an input sequence
$x_{in}=\\{x_{0},x_{1},...,x_{\tau}\\}$ is fed in to produce a sequence of
outputs
$x^{\prime}_{out}=\\{x^{\prime}_{\tau+1},x^{\prime}_{\tau+2},...,x^{\prime}_{\tau+k}\\}$.
Prediction losses are calculated with respect to some metric $d$ and the real
output $x_{out}=\\{x_{\tau+1},x_{\tau+2},...,x_{\tau+k}\\}$,

$\mathcal{L}=d(x^{\prime}_{out},x_{out}).$ (4)

The gradient of the loss taken with respect to node weights is calculated by
propagating the error backward through time steps and used to update node
weights (see. Figure 1).

Figure 1: Schematic of typical backpropagation through time calculation for
training recurrent neural networks with an input length of 4. Input time
series $x_{t}$ is used to evolve node states $s_{t}$ during the forward pass.
Node states are used to calculate a value for the error $E_{t}$ at a given
time $t$. Derivatives are taken with respect to node states and weights.
Variables $c_{n}$ correspond to the weight of the link joining $s_{n-1}$ to
$s_{n}$ and is contained within $C^{h}(t)$

.

From observation, the hidden layer of an RNN functions identically to a
dynamical network. We aim to investigate if the BPTT approach can be adapted
to regress the weights of a dynamical network given a set of input node
states.

##  3 Methods

This paper proposes a novel method of simultaneously inferring the local
dynamics $f$ and connectivity structure $C$ of a dynamical network with a
known coupling function $g$ by only observing node states. The proposed method
consists of three main stages, namely initialisation, backpropagation and
decoupling. A flowchart outlining this process is provided (see Figure 2).

Figure 2: Overview of the steps in the proposed backpropagation regression
method. Components being trained (red) in each step. Regressed components are
given with solid lines.

###  3.1 Backpropagation Regression Algorithm

#### Initialisation

The initialisation stage consists of constructing a local model $\hat{f}$ that
approximates the true local dynamics $f$. Due to the complex interplay between
the local $f$ and coupling $g$ dynamics, getting an accurate coupling-free
estimation of $f$ directly from node signals is not possible.

The estimated local model $\hat{f}$ is defined as a map from state space
$\mathbf{x}(t)$ to vector field $\mathbf{\dot{x}}(t)$. The time series is
numerically differentiated in order to retrieve $\mathbf{\dot{x}}(t)$,

$\mathbf{\dot{x}}(t)\approx\frac{\mathbf{\dot{x}}(t+\delta
t)-\mathbf{\dot{x}}(t)}{\delta t}.$ (5)

For this stage of the method, we assume that the coupling effects from $g$ is
relatively small compared to local dynamics $f$. This allows a mean field
approach to be taken when constructing an approximate local model $\hat{f}$.
We note that this method requires that the ‘attractor’ (coupling with other
neighbours in the dynamical network may cause spurious intersections of
trajectories) is not too distorted from the attractor of the non-coupled case.
Here, we calculate an estimate of the true vector field at each point in state
space as an average over its $K$ nearest neighbours (see Figure 3),

$\mathbf{\hat{\dot{x}}}(t)=\frac{1}{K}\sum_{i=1}^{K}\frac{\mathbf{\dot{x}}^{(i)}(t+\delta
t)-\mathbf{\dot{x}}^{(i)}(t)}{\delta t},$ (6)

where $\mathbf{\dot{x}^{(i)}}(t)$ corresponds to the velocity of the $i^{th}$
closest neighbour. The identification of neighbours are done with respect to
the Euclidean distance with all points in the time series (i.e at all times
$t$) being considered potential candidates of nearest neighbours.

Figure 3: Mean field approach where effective vector field (blue) at a point
in phase space is taken as the average of $K$ nearest neighbours’ vector field
influenced with neighbouring node interactions.

The estimated local dynamics $\hat{f}$ is defined as a mapping
$\hat{f}:\mathbf{x}(t)\to\mathbf{\hat{\dot{x}}}(t)$. A feedforward network
with 2 hidden layers was used to learn the mapping due to its simplicity and
ease of implementation. Each hidden layer is composed of 128 neurons with a
sigmoid activation function. This estimated local dynamics $\hat{f}$
calculated from the mean field approach acts as an initial proxy for the true
local dynamics, which will then be further refined to remove coupling effects
on the vector field during the backpropagation and retraining stages of the
algorithm. In addition to the local dynamics $\hat{f}$, values for the
estimated coupling adjacency matrix $\hat{C}$ must also be initialised.
However, there is no restriction on what values of $\hat{C}$ may be taken.

#### Backpropagation

The main aim of the backpropagation stage is to provide an improved estimate
for the coupling adjacency matrix $\hat{C}$. Similar to the BPTT algorithm, a
forward pass consisting of freerun predictions is first calculated using the
estimated $\hat{f}$ and $\hat{C}$ with randomly selected initial values from
the observed time series in order to properly sample across the whole
attractor,

$\displaystyle\hat{\dot{\mathbf{x}}}_{i}=\hat{f}(\mathbf{x}_{i})+\sum_{i\neq
j}\hat{c}_{i,j}g(\hat{\mathbf{x}}_{i},\hat{\mathbf{x}}_{j}),$ (7a)
$\displaystyle\hat{{\mathbf{x}}}_{i}(t+\delta
t)=\hat{{\mathbf{x}}}_{i}(t)+\delta t\;\hat{\dot{\mathbf{x}}}_{i}(t).$ (7b)

For simplicity, we assume coupling only in the first component of each time
series.

Freerun predictions for all nodes $\mathbf{\hat{x}_{i}}(t)$ are calculated by
recursively evaluating Equations 7. The loss in freerun predictions due to
errors in the local model $\hat{f}$ and network weights $\hat{C}$ can be
calculated using some distance function. Here, $n$ is the length of the
freerun prediction, $N$ is the number of nodes. This loss at each node,
$E_{i}$, is calculated using the $L_{2}$-norm across the total freerun
trajectory {$\hat{\mathbf{x}}_{i}(t_{0}),...,\hat{\mathbf{x}}_{i}(t_{n})$},

$\mathcal{L}=\sum_{t=t_{0}}^{t=t_{n}}\sum_{i=1}^{N}E_{i}(t)=\sum_{t=t_{0}}^{t=t_{n}}\sum_{i=1}^{N}\lVert\hat{\mathbf{x}}_{i}(t)-\mathbf{x}_{i}(t)\rVert^{2}.$
(8)

For regressing network weights, the loss gradient with respect to $\hat{C}$
must be calculated and used to update the estimate of $\hat{C}$,

$\medmath{\frac{d\mathcal{L}}{d\hat{C}}=\sum_{t=t_{0}}^{t=t_{n}}\sum_{i=1}^{N}\frac{\partial
E_{i}(t)}{\partial\hat{C}}=\sum_{t=t_{0}}^{t=t_{n}}\sum_{i=1}^{N}\sum_{d}2(\hat{\mathbf{x}}^{(d)}_{i}(t)-\mathbf{x}^{(d)}_{i}(t))\frac{\partial\hat{\mathbf{x}}^{(d)}_{i}(t)}{\partial\hat{C}}.}$
(9)

This requires the estimation of partial derivatives with respect to the
coupling weights $c_{ij}$. The recursive expressions for these derivatives are
given in the Appendix. To provide more stability in the regression, the
gradient for each step of backpropagation is averaged over the multiple
randomly chosen initial inputs $\mathbf{x}(t_{0})$. An overview of the
modified backpropagation algorithm is given in Figure 4.

Figure 4: Schematic of the two node modified backpropagation algorithm for
calculating weight gradients. Forward pass (blue and black), and
backpropagation (red).

#### Decoupling and Retraining

The use of a mean field approximation to construct an initial local model
$\hat{f}$ results in significant errors in the estimation of the vector field
for the coupled dimensions. These errors limit the ability for the algorithm
to estimate the true coupling weights as accurate regressions assume that
$f=\hat{f}$. By using estimated coupling values $\hat{C}$ to partially
decouple calculated values of the vector field, the output space of $\hat{f}$
changes to better reflect the true local dynamics $f$ (i.e. without coupling
effects). To adapt to this change, the model must be retrained on data of the
new partially decoupled output space. The process of decoupling the vector
field and truncating coupling interaction effects allows for gradual
correction of the errors introduced from the mean field approximation of the
vector field.

After attaining an improved estimate of the weighted adjacency matrix
$\hat{C}$, the calculated coupling weights can be used to partially decouple
the observed node time series according to Equation 10.

$\hat{f}(x_{i}(t_{n}))\approx\frac{x_{i}(t_{n+1})-x_{i}(t_{n})-\delta
t\;\sum_{i\neq j}\hat{c}_{ij}g(x_{i}(t_{n}),x_{j}(t_{n}))}{\delta t}.$ (10)

The new calculated value of the vector field $\hat{f}$ can then be used to
improve the previously constructed local dynamics model. The final decoupling
step is alternated with the backpropagation step to gradually improve both the
local model and estimate weighted adjacency matrix $\hat{C}_{ij}$. A much
lower model learning rate $\eta^{\prime}<\eta$ must be used to ensure that the
structure of the initial model is preserved. Using a small learning rate also
has the added benefit of providing additional stability to the algorithm as it
minimises drastic changes in the loss landscape that can affect convergence.

###  3.2 Performance Assessment

We need to introduce several measures to assess performance for the proposed
backpropagation method. Firstly, prediction horizons were measured for the
regressed dynamical network to assess the overall quality of the combined
model inclusive of local dynamics and coupling weights. The prediction horizon
for node $i$ was defined as the time $t_{p,i}$ where prediction errors
exceeded a set threshold $\epsilon$. The overall prediction horizon $t_{p}$
for the entire dynamical network is taken as the average prediction horizon
across all nodes $i$,

$t_{p}=\langle
t_{p,i}\;|\;\lVert\hat{\mathbf{x}}_{i}(t_{p,i})-\mathbf{x}_{i}(t_{p,i})\rVert^{2}<\epsilon\rangle_{i=1}^{N}.$
(11)

The norm of the weight error matrix, $\epsilon_{C}$ is used to assess the
overall error in the regressed weighted adjacency matrix and is normalised the
by the true weight values,

$\epsilon_{C}=\frac{\lVert\hat{C}-C\rVert^{2}}{\lVert C\rVert^{2}}.$ (12)

Mutual information $I(X,Y)$ is a nonlinear measure of the correlation between
two random variables $X$ and $Y$. In simple terms, describes the amount of
information that is shared between any two random variables $X$ and $Y$. Low
mutual information implies independence between $X$ and $Y$. Mutual
information is usually calculated as a function of marginal and joint
probabilities,

$I(X,Y)=\int_{\mathcal{X}}\int_{\mathcal{Y}}p_{(X,Y)}(x,y)\ln\left({\frac{p_{(X,Y)}(x,y)}{p_{X}(x)p_{Y}(y)}}\right)\,dydx.$
(13)

The quality of a local model can be accessed by comparing the mutual
information between the true states of a trajectory $\mathbf{x}(t)$, and those
predicted by the model $\hat{\mathbf{x}}(t)$ at each point in time $t$. Hence,
we can define the average mutual information $I(\tau)$ with respect to some
chosen freerun prediction length $\tau$,

$I(\tau)=\frac{1}{\tau}\sum_{t=1}^{\tau}p(\mathbf{x}(t),\hat{\mathbf{x}}(t))\ln\left({\frac{p(\mathbf{x}(t),\hat{\mathbf{x}}(t))}{p(\mathbf{x}(t))p(\hat{\mathbf{x}}(t))}}\right),$
(14)

where $p(\mathbf{x}(t)),p(\hat{\mathbf{x}}(t))$ and
$p(\mathbf{x}(t),\hat{\mathbf{x}}(t))$ are the joint and marginal
probabilities of observing (or predicting for the case of $\mathbf{\hat{x}}$)
a given state at time $t$. For each value of $\tau$ in $I(\tau)$, the
distributions of $p(\mathbf{x}(t)),p(\hat{\mathbf{x}}(t))$ and
$p(\mathbf{x}(t),\hat{\mathbf{x}}(t))$ are calculated empirically using
histograms based on observations of $\mathbf{x}(t)$ and $\mathbf{\hat{x}}(t)$
for $t\in[1,\tau]$.

The function $I(\tau)$ can be interpreted as measure of the prediction
performance of the local model for a prediction length of $\tau$ steps. For
chaotic systems, $I(\tau)$ of an imperfect model will be a decreasing
function. The contour of $I(\tau)$ tracks the collapse of the local model’s
prediction performance for increasing prediction lengths $\tau$.

The performance of a trained local model $\hat{f}$ can be assessed by using it
to produce freerun predictions of separate system with only a single node
(i.e. no external interactions). A score $S$ related to the mutual information
between the trajectories from the real local model $\mathbf{x}(t)$ and trained
local model $\hat{\mathbf{x}}(t))$ can be calculated using the cumulative area
under the average mutual information curve,

$S=\int_{1}^{\tau_{max}}I(\tau)d\tau,$ (15)

where larger values of $S$ indicate more persistent mutual information and a
better predictive local model.

The mutual information based score $S$ is only calculated with respect to only
one component of the predicted time series (i.e. $x$ component). This is
sufficient for the systems analysed as all components of the time series
($x,y,z$) feedback to each other. Hence, collapse in the prediction
performance of one component will quickly result in poor predictions in all
other components. Additionally, $I(\tau)$ for each component cannot be summed
together as it assumes that all components evolve independently.

To assess the mutual information of the combined model of both regressed
weights and trained local model, $S$ is calculated as the sum of the scores
across all nodes. Whilst nodes are not independent, the sparsity of network
connections mean that collapse in the prediction performance of one node does
necessarily propagate quickly to all other nodes in the network.

##  4 Results

Three tests were done to assess the performance of the proposed adapted
backpropagation regression algorithm. These were the cases of the small 16
node network, large 64 node network, and networks with non-positive asymmetric
coupling weights. All tests were conducted on random networks with each edge
existing with probability $p=log(N)/N$ to minimise the likelihood of disjoint
networks. Edges are bidirectional in all networks except the asymmetric
network case. Values for the node coupling strengths were normally distributed
$(\mu,\sigma)=(0.15,0.02)$. The mean coupling strength $\mu$ is assumed to be
non-zero to ensure that regressed values clearly distinguish between the
absence and presence of edges.

Input data of $25000$ steps with an additional $2000$ steps washout were
generated with timestep $dt=0.02$ using the standard RK4 algorithm. Washout
steps were discarded exclude any transient dynamics. An initial learning rate
of $\eta=0.001$ was used for training the feedforward model during
initialisation. A lower learning rate of $\eta^{\prime}=0.0002$ was used for
all subsequent retraining iterations with each training period consisting of
$30$ epochs. A preset number of $40$ or $80$ refit iterations were chosen
before ending regression. However, we note that it may be possible to
implement a stopping criterion based on plateauing of either the prediction
errors or regressed weight matrix $\hat{C}$.

Performance measures $(S,t_{p})$ for each model were evaluated over 8 randomly
chosen initial conditions.To assess the performance of regressed models, each
performance measure was compared against a control case where the entire
network structure and local model is fully known but initial conditions are
perturbed by $\xi\sim N(0,0.005^{2})$ where $\xi$ is in scaled units. This
method was chosen over the calculation of dynamical invariants such as
Lyapunov exponents due to its ease of computation, especially in high
dimensional systems such as dynamical networks. Due the chaotic nature of the
system’s oscillators, small deviations in initial conditions will grow
exponentially. Deviations the initial conditions that are sufficiently small
are able to track nearby trajectories for some period of time before
diverging. Therefore, scores that match those of the control case imply that
the regressed models are able to perform almost as well as the case where the
system dynamics are fully known, but contaminated with noise in its initial
conditions.

###  4.1 Positive Symmetric Networks

The fully observable case assumes that all node time series are accessible for
regression. Gradient losses $\mathcal{L}$ are calculated at each
backpropagation step and used to regress the network adjacency matrix
$\hat{C}$. The regression algorithm was run for a total of $30$ refit
iterations each consisting of $300$ backpropagation steps followed by
decoupling and retraining of the model. A total of 8 randomly generated 16
node networks were tested. The quality of the regressed weights and local
model for each network was assessed by running freerun predictions with 16
randomly generated initial conditions.

Figure 5: Regression results for the Lorenz 16 node system normally
distributed coupling weights $(\mu,\sigma^{2})=(0.15,0.02)$ and connection
probability $p$ over 30 iterations. First two data points correspond to the
values attained during initialisation. Figure 6: Regression of Lorenz 16 node
network weights (top) at different number of iterations compared to the
normalised error in the weights (bottom). True weights (right) are given for
comparison showing good agreement with the regressed results. Figure 7:
Regression results for the Chua 16 node system operating in the single scroll
regime. Coupling weights are normally distributed with
$(\mu,\sigma^{2})=(0.15,0.02)$ and connection probability $p$ over 30
iterations. First two data points correspond to the values attained during
initialisation. Figure 8: Regression of Chua 16 node network weights (top) at
different number of iterations compared to the normalised error in the weights
(bottom). True weights (right) are given for comparison showing good agreement
with the regressed results.

The backpropagation method was able to accurately regress the weights of both
the Lorenz (Figure 5) and Chua (Figure 7) dynamical networks. The average
normalised error in the regressed weights for these systems at the end of 30
iterations was 9% and 7% for the Lorenz and Chua networks respectively. The
constructed trained local model for both systems also showed a large
improvement in $S$ with 71% increase for Lorenz and a 38% increase for Chua
when compared to the initial local model based purely on the mean field
approach (i.e. Iteration 0). The mutual information of the combined local
model and coupling structure increased by 65% and 59% for the Lorenz and Chua
networks respectively. The prediction horizons for both systems were found to
increase with each additional iteration. Both Chua and Lorenz networks were
found to perform as well or better than the control case where $\xi=0.005$.

###  4.2 Large 64-node Network

A single larger system of 64 random coupled Lorenz oscillators was used to
test the scalability of the system. The algorithm was run for 80 iterations to
account for slightly slower convergence. Performance metrics identical to the
regular 16 node case was used with the network tested against 16 different
initial conditions (see Figures 9 and 10).

We find that the backpropagation method was also able to accurately regress
the desired weights of the the coupled system with a final normalised weight
error of 9%. The mutual information of the local model and combined network
system was also found to increase by 51% and 81% by the end of the iterations
respectively.

Figure 9: Regression results for the larger Lorenz 64 node system normally
distributed coupling weights $(\mu,\sigma^{2})=(0.15,0.02)$ and connection
probability $p$ over 80 iterations. First two data points correspond to the
values attained during initialisation. Figure 10: Regression of Lorenz 64 node
network weights (top) at different number of iterations compared to the
normalised error in the weights (bottom). True weights (right) are given for
comparison showing good agreement with the regressed results.

###  4.3 Negative Asymmetric Coupling

Many dynamical systems assume the case of positive coupling resulting in the
attraction of nearby trajectories. In addition to this common case, we
consider dynamical systems with both negative and positive coupling where
repulsive behaviour between nodes is possible. Additionally, we also consider
the effects where the adjacency matrix is not symmetrical. To do this, a 16
node adjacency matrix was generated with each non diagonal entry being non-
zero with a probability $p=log(N)/N$ and is not guaranteed to be symmetric.
However, each node has a 0.25 probability of being changed to a negative
value. The time series was then generated similarly to the original test case
with the Lorenz equations (see Figures 11 and 12).

Figure 11: Regression results for the Lorenz 16 node system normally
distributed coupling weights $(\mu,\sigma^{2})=(0.15,0.02)$ and connection
probability $p$ over 80 iterations. Weights are asymmetric and changed to its
negative value with probability of 0.25. First two data points correspond to
the values attained during initialisation. Figure 12: Regression of Lorenz 16
node network with negative and positive weights (top) at different number of
iterations compared to the normalised error in the weights (bottom). True
weights (right) are given for comparison showing good agreement with the
regressed results.

The backpropagation method was found to also work equally well in regressing
cases with negative weights and asymmetric coupling. The results for the
network prediction horizon, log weight errors and mutual information measures
show a similar trend to their symmetric, positively weighted counterparts with
a final normalised weight error of 5%, and an 58% relative increase in mutual
information for the constructed local model. The mutual information of the
combined network was found to increase by 75% when compared to the first
iteration.

###  4.4 Noise Effects

In order to test the robustness of the method, the backpropagation approach
was applied to time series from a 16 node Lorenz dynamical network with
additive Gaussian noise to simulate measurement noise. The Gaussian noise term
$\epsilon_{\xi}(t)\sim N(0,\xi)$ was added after normalising the data to zero
mean and unit standard deviation. Because backpropagation regression requires
the evaluation of numerical derivatives, the noisy input data was first
smoothed using spline regression [40, 41] with a regularisation parameter
$\lambda=10$. The usage of cubic splines also has the advantage or producing
relatively smooth derivatives.

The noise magnitude was varied with $\xi\in[0,0.1]$ with the maximum $\xi$
value corresponding to a signal-noise ratio of 20dB. Backpropagation
regression was run for 40 iterations in line with those for the noiseless
Lorenz and Chua 16 node systems. Each configuration was repeated for 8
randomly chosen network and initial conditions.

The effect of noise hampers performance and results in the addition of
spurious incorrect weights during regression (see Figure 13). To address this,
a threshold value of 0.04 is selected to filter the weights following
regression. Regressed weights that do not exceed this value are set to 0. The
threshold value is visually selected from the histogram distribution of
weights. As the network weights are normally distributed with non-zero mean,
$c_{ij}\sim N(0.15,0.02)$, true weights will tend to have regressed value much
greater than zero. In contrast, noise effects will yield to low valued
spurious weights. This results in a bimodal distribution of weights. Hence,
the threshold value can be selected as a value that distinguishes between the
two modes.

The performance of the regression algorithm weights was found to be relatively
robust against noise. For low noise levels $\xi\leq 0.04$, the effects of
noise appear follow a linear relationship (see Figure 13). The increase in
weight errors $\epsilon_{C}$ slow for higher levels of noise. In the latter
case, the weight filtration was found to be very effective in removing
spurious edge weights allowing a further decrease in weight errors. However,
we found that higher values of noise with $\xi\geq 0.15$ resulted in
instability of the regression algorithm and divergence in the weight errors.

The effect of noise was found to significantly affect the quality of the local
dynamics models (see Figure 14). For small noise levels $\xi\leq 0.01$, the
approximations $\hat{f}$ of the local dynamics was relatively stable. However,
there was a large difference in prediction performance between the noiseless
and noisy case. Noise effects are more detrimental for $\xi>0.01$ where the
long run predictions of the resulting models no longer conform to the
attractor of the original Lorenz system. In several instances, it was found
that the models slowly drifted and converged to a fixed point. For $\xi\geq
0.07$, the observational noise results in a local model that is worse than the
initial mean field estimate.

The poor learning of the local model can be attributed to the combined effects
of noise changing the temporal structure of the time series and the relatively
low sampling rates of the time series. An integration step of $dt=0.02$
corresponds to approximately 35 points per oscillation. Each oscillation
conists of a short transient extrema. Additive noise result in over or
underestimation of these peaks, which can have a detrimental effect on
calculated vector field. The noise tests were repeated at a lower integration
step of $dt=0.002$, which was smoothed and downsampled to $dt=0.02$. This
configuration allowed a better treatment of the additive noise and yielded
more robust results.

Overall, we find that the performance of this method under noise is comparable
to other methods used to analysed dynamical networks. For local dynamics, the
attractor of the regressed local model dynamics began to collapse at
$\xi=0.05$ corresponding to a signal-noise ratio of 26dB. This is comparable
to the test results conducted by [31] where the methods SINDy, ARNI and Gao’s
method exhibit large decreases in performance at approximately 55dB, 40dB and
25dB respectively.

Figure 13: Weight errors $\epsilon_{C}$ after 40 iterations for
backpropagation regression applied to 16 node Lorenz dynamical networks with
varying levels of additive Gaussian noise. Each noise configuration was
repeated with 8 randomly chosen networks and initial conditions. Error bars
correspond to interquartile range of weight errors. Weights with a magnitude
less than a threshold 0.04 are truncated and set to zero. The threshold value
selected selected from the histogram distribution of all weights. Figure 14:
Mutual information of the local models for configurations with increasing
observational noise $\xi$. Model scores are compared at the beginning (blue)
and end (red) of the 40 refit iterations. Mutual information is compared
against the control case where the exact model known with perturbed initial
conditions.

##  5 Method Limitations

There are several data and system requirements needed for the application of
the backpropagation regression algorithm. Firstly, the backpropagation
regression method requires networks that are not fully synchronised. This is
because differences in the trajectories of coupled nodes are used to decouple
observed node signals to infer links. Networks operating in the fully
synchronised regime will have fully coincident trajectories and cannot be used
to infer weights.

Secondly, reliance of a mean field approximation during the initialisation
stage assumes that node trajectories do not deviate too much from the
attractor of the local dynamics of the uncoupled single node case. Small
deviations will inevitably occur, and is accounted for during the decoupling
and retraining stage of each iteration. However, deviations that are too large
may cause convergence problems.

Thirdly, there must be sufficient full node state observations with
trajectories that fully explore the state space for the regression to converge
and produce a good local model. The method was found to be convergent with
time series lengths as low as 5000 time steps with $dt=0.02$, excluding an
initial washout period of $2000$ time steps to account for dynamical
transients.

##  6 Conclusion

Dynamical networks are an interesting and useful framework to analyse complex
systems with multiple interacting components. However, application of this
framework is usually difficult due to lack of information pertaining to
features such as the connectivity structure and local dynamics. As a result,
there is a need for a method to extract these features purely from observed
node states.

In this paper, we propose an adapted backpropagation method for inferring the
local model and connectivity structure of dynamical networks purely from node
time series. This method draws inspiration from the backpropagation through
time (BPTT) algorithm commonly used to train RNNs and adapts it for the
purpose of regressing couplings weights in a dynamical network. Regressed
values for the network coupling are then used to decouple the observed node
time series and improve the construction of the local model without coupling
effects. The two steps of backpropagation weight regression and decoupling
alternate until sufficient convergence is achieved yielding a final local
model and regressed coupling weights.

The method was tested on dynamical networks of nodes of identical chaotic
oscillators with successful results. Two types of chaotic oscillators, Lorenz
and single-scrolled Chua, were tested. The backpropagation regression method
was able to accurately identify the coupling weights down to a margin of 7%.
The regressed models were also compared to a control case where the true local
model is known, but whose initial states are perturbed by a small error. It
was found that the trained local model was comparable to the control case.
Combined, the regressed weights and constructed local model were able to
achieve respectable prediction horizons for the 16 node dynamical network
compared to a control case. Despite its poor scalability, the method was also
able to accurately reconstruct the local model and coupling weights of a
larger 64 node network with similar results to the smaller 16 node
counterparts. Additional testing was also done on dynamical networks with
negative and asymmetric coupling with similar performance.

The robustness of the method was also tested against noise and network
heterogeneity (see Appendix). The algorithm was found to be convergent for
small to moderate amounts of observational noise. However, larger noise levels
resulted in poor local models and potential instability. In contrast, the
presence of network heterogeneity did not significantly affect the performance
of the backpropagation regression algorithm.

We acknowledge that the method proposed in this paper is a proof of concept
that, whilst functional, has multiple areas of further research and potential
refinement. Firstly, there remains further work on the refinement of the speed
and complexity of the algorithm to be more scalable. Secondly, there remains a
need for further testing the on the effect of hyperparameters and
initialisation on algorithm performance. Thirdly, it would be beneficial to
investigate the performance of this method on different dynamical network
structures. There also remains the possibility of investigating how other
machine learning methods used in recurrent neural networks may be applied to
the backpropagation regression algorithm.

## Acknowledgments

E.T. is supported by a Robert and Maude Gledden Postgraduate Research
Scholarship and Australian Government RTP Scholarship at The University of
Western Australia. M.S. and D.C.C. acknowledge the support of the Australian
Research Council through the Centre for Transforming Maintenance through Data
Science (grant number IC180100030), funded by the Australian Government.

## References

  * [1] D. M. Abrams and S. H. Strogatz, “Chimera states for coupled oscillators,” _Physical Review Letters_ , vol. 93, no. 17, p. 174102, 2004.
  * [2] A. Andreev, N. Frolov, A. Pisarchik, and A. Hramov, “Chimera state in complex networks of bistable hodgkin-huxley neurons,” _Physical Review E_ , vol. 100, no. 2, p. 022224, 2019.
  * [3] Y. Zhu, Z. Zheng, and J. Yang, “Chimera states on complex networks,” _Physical Review E_ , vol. 89, no. 2, p. 022914, 2014.
  * [4] S. Majhi, B. K. Bera, D. Ghosh, and M. Perc, “Chimera states in neuronal networks: a review,” _Physics of Life Reviews_ , vol. 28, pp. 100–121, 2019\.
  * [5] E. M. Izhikevich, _Dynamical Systems in Neuroscience_. MIT press, 2007.
  * [6] B. Wu, D. Zhou, F. Fu, Q. Luo, L. Wang, and A. Traulsen, “Evolution of cooperation on stochastic dynamical networks,” _PLoS One_ , vol. 5, no. 6, p. e11187, 2010.
  * [7] T. Liu, J. Zhao, and D. J. Hill, “Exponential synchronization of complex delayed dynamical networks with switching topology,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 57, no. 11, pp. 2967–2980, 2010\.
  * [8] A. R. Hota, T. Sneh, and K. Gupta, “Impacts of game-theoretic activation on epidemic spread over dynamical networks,” _SIAM Journal on Control and Optimization_ , no. 0, pp. S92–S118, 2021.
  * [9] M. Rosenblum, M. Frühwirth, M. Moser, and A. Pikovsky, “Dynamical disentanglement in an analysis of oscillatory systems: an application to respiratory sinus arrhythmia,” _Philosophical Transactions of the Royal Society A_ , vol. 377, no. 2160, p. 20190045, 2019.
  * [10] S. G. Shandilya and M. Timme, “Inferring network topology from complex dynamics,” _New Journal of Physics_ , vol. 13, no. 1, p. 013004, 2011.
  * [11] G. Stepaniants, B. W. Brunton, and J. N. Kutz, “Inferring causal networks of dynamical systems through transient dynamics and perturbation,” _Physical Review E_ , vol. 102, no. 4, p. 042309, 2020.
  * [12] A. Banerjee, J. Pathak, R. Roy, J. G. Restrepo, and E. Ott, “Using machine learning to assess short term causal dependence and infer network links,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 29, no. 12, p. 121104, 2019.
  * [13] A. Banerjee, J. D. Hart, R. Roy, and E. Ott, “Machine learning link inference of noisy delay-coupled networks with optoelectronic experimental tests,” _Physical Review X_ , vol. 11, no. 3, p. 031014, 2021.
  * [14] D. Napoletani and T. D. Sauer, “Reconstructing the topology of sparsely connected dynamical networks,” _Physical Review E_ , vol. 77, no. 2, p. 026103, 2008.
  * [15] C. Weistuch, L. Agozzino, L. R. Mujica-Parodi, and K. A. Dill, “Inferring a network from dynamical signals at its nodes,” _PLoS Computational Biology_ , vol. 16, no. 11, 2020.
  * [16] M. Paluš, “Coupling in complex systems as information transfer across time scales,” _Philosophical Transactions of the Royal Society A_ , vol. 377, no. 2160, p. 20190094, 2019.
  * [17] M. Kornilov, I. Sysoev, D. Astakhova, D. Kulminsky, B. Bezruchko, and V. Ponomarenko, “Reconstruction of the coupling architecture in the ensembles of radio-engineering oscillators by their signals using the methods of granger causality and partial directed coherence,” _Radiophysics and Quantum Electronics_ , vol. 63, no. 7, pp. 542–556, 2020.
  * [18] X. Wu, W. Wang, and W. X. Zheng, “Inferring topologies of complex networks with hidden variables,” _Physical Review E_ , vol. 86, no. 4, p. 046106, 2012\.
  * [19] F. Takens, “Detecting strange attractors in turbulence,” in _Dynamical Systems and Turbulence, Warwick 1980_. Springer, 1981, pp. 366–381.
  * [20] J.-S. Zhang and X.-C. Xiao, “Predicting chaotic time series using recurrent neural network,” _Chinese Physics Letters_ , vol. 17, no. 2, p. 88, 2000\.
  * [21] M. Small and C. K. Tse, “Minimum description length neural networks for time series prediction,” _Physical Review E_ , vol. 66, no. 6, p. 066701, 2002\.
  * [22] H. Jaeger, “The “echo state” approach to analysing and training recurrent neural networks-with an erratum note,” _Bonn, Germany: German National Research Center for Information Technology Gesellschaft für Mathematik und Datenverarbeitung mbH (GMD) Technical Report_ , vol. 148, p. 13, 2001.
  * [23] J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, “Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 27, no. 12, p. 121102, 2017.
  * [24] Z. Lu, B. R. Hunt, and E. Ott, “Attractor reconstruction by machine learning,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 28, no. 6, p. 061104, 2018.
  * [25] A. Haluszczynski and C. Räth, “Good and bad predictions: Assessing and improving the replication of chaotic attractors by means of reservoir computing,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 29, no. 10, p. 103143, 2019.
  * [26] E. Tan, D. Corrêa, T. Stemler, and M. Small, “Grading your models: Assessing dynamics learning of models using persistent homology,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 31, no. 12, p. 123109, 2021.
  * [27] M. J. Panaggio, M.-V. Ciocanel, L. Lazarus, C. M. Topaz, and B. Xu, “Model reconstruction from temporal data for coupled oscillator networks,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 29, no. 10, p. 103116, 2019.
  * [28] W.-X. Wang, Y.-C. Lai, and C. Grebogi, “Data based identification and prediction of nonlinear and complex dynamical systems,” _Physics Reports_ , vol. 644, pp. 1–76, 2016.
  * [29] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016.
  * [30] W.-X. Wang, R. Yang, Y.-C. Lai, V. Kovanis, and M. A. F. Harrison, “Time-series–based prediction of complex oscillator networks via compressive sensing,” _EPL (Europhysics Letters)_ , vol. 94, no. 4, p. 48006, 2011.
  * [31] T.-T. Gao and G. Yan, “Autonomous inference of complex network dynamics from incomplete and noisy data,” _Nature Computational Science_ , vol. 2, no. 3, pp. 160–168, 2022.
  * [32] J. Casadiego, M. Nitzan, S. Hallerberg, and M. Timme, “Model-free inference of direct network interactions from nonlinear collective dynamics,” _Nature Communications_ , vol. 8, no. 1, pp. 1–10, 2017.
  * [33] D. Eroglu, M. Tanzi, S. van Strien, and T. Pereira, “Revealing dynamics, communities, and criticality from data,” _Physical Review X_ , vol. 10, no. 2, p. 021047, 2020.
  * [34] P. J. Werbos, “Backpropagation through time: what it does and how to do it,” _Proceedings of the IEEE_ , vol. 78, no. 10, pp. 1550–1560, 1990.
  * [35] E. N. Lorenz, “Deterministic nonperiodic flow,” _Journal of Atmospheric Sciences_ , vol. 20, no. 2, pp. 130–141, 1963.
  * [36] J. Kengne, “On the dynamics of chua’s oscillator with a smooth cubic nonlinearity: Occurrence of multiple attractors,” _Nonlinear Dynamics_ , vol. 87, no. 1, pp. 363–375, 2017.
  * [37] S. Raubitzek and T. Neubauer, “Taming the chaos in neural network time series predictions,” _Entropy_ , vol. 23, no. 11, p. 1424, 2021.
  * [38] P. Dubois, T. Gomez, L. Planckaert, and L. Perret, “Data-driven predictions of the lorenz system,” _Physica D: Nonlinear Phenomena_ , vol. 408, p. 132495, 2020.
  * [39] M. Sangiorgio and F. Dercole, “Robustness of LSTM neural networks for multi-step forecasting of chaotic time series,” _Chaos, Solitons and Fractals_ , vol. 139, p. 110045, 2020.
  * [40] P. J. Green and B. W. Silverman, _Nonparametric regression and generalized linear models: a roughness penalty approach_. CRC Press, 1993.
  * [41] C. H. Reinsch, “Smoothing by spline functions,” _Numerische Mathematik_ , vol. 10, no. 3, pp. 177–183, 1967.

Eugene Tan is a Robert and Maude Gledden scholar and PhD candidate at the
University of Western Australia. He completed his B.Sc and MPE in mechanical
engineering at the University of Western Australia.
---
Dèbora Corrêa completed her PhD degree in Applied Physics with a focus on
Computational Physics at the University of Sao Paulo (USP). She is currently a
lecturer at the University of Western Australia
---
Thomas Stemler completed his Ph.D. at the Technische Universität Darmstadt. He
is currently a senior lecturer at the University of Western Australia.
---
Michael Small is a professor and CSIRO-UWA Chair in Complex Engineering
Systems at the University of Western Australia. He completed his Ph.D. in
applied mathematics focusing on nonlinear dynamical systems and time series
analysis at The University of Western Australia.
---