# An Analysis Of Protected Health Information Leakage
In Deep-Learning Based De-Identification Algorithms

Salman Seyedi, 1 Li Xiong, 2 Shamim Nemati, 3 Gari D. Clifford 1,4

###### Abstract

The increasing complexity of algorithms for analyzing medical data, including
de-identification tasks, raises the possibility that complex algorithms are
learning not just the general representation of the problem, but specifics of
given individuals within the data. Modern legal frameworks specifically
prohibit the intentional or accidental distribution of patient data, but have
not addressed this potential avenue for leakage of such protected health
information. Modern deep learning algorithms have the highest potential of
such leakage due to complexity of the models. Recent research in the field has
highlighted such issues in non-medical data, but all analysis is likely to be
data and algorithm specific. We, therefore, chose to analyze a state-of-the-
art free-text de-identification algorithm based on LSTM (Long Short-Term
Memory) and its potential in encoding any individual in the training set.
Using the i2b2 Challenge Data, we trained, then analyzed the model to assess
whether the output of the LSTM, before the compression layer of the
classifier, could be used to estimate the membership of the training data.
Furthermore, we used different attacks including membership inference attack
method to attack the model. Results indicate that the attacks could not
identify whether members of the training data were distinguishable from non-
members based on the model output. This indicates that the model does not
provide any strong evidence into the identification of the individuals in the
training data set and there is not yet empirical evidence it is unsafe to
distribute the model for general use.

An electronic health record, or electronic medical record (EMR) includes a
wealth of information in the form of both physiological data and structured or
free text. The latter is often replete with protected health information (PHI)
and personal identifiable information (PII). As a result, there has been much
attention paid to the notion of secure processing and sharing. The integrity
of patients' personal information and related privacy are governed in the US
by the Health Insurance Portability and Accountability Act of 1996 (HIPAA).
Although intended to make medical data portable, it has had a much larger
effect in ensuring privacy through de-identification of shared data. The HIPAA
privacy rules provide two avenues that one must follow to meet the de-
identification standards: ‘expert determination’ and ‘safe harbor’. The safe
harbor de-identification method requires the removal of Protected Health
Information (PHI) which is a list of 18 categories of identifiers. The process
of de-identification of a particular document can be performed using different
means, but based on the enormous amount of EMR which grows every moment, the
sheer volume not only greatly incentivizes the automation of the process as
much as possible, but one can argue any pragmatic approach have to rely
intensively on utilizing power of computers. Because of this, the automation
of the de-identification process has been of great interest and different
algorithms have been developed over the years to facilitate the de-
identification by automating the finding/labeling of the PHI (Neamatullah et
al. 2008; Yang and Garibaldi 2015; Dernoncourt et al. 2016).

In recent years due to improvements in the deep learning architectures and
recurrent neural networks (RNN), not only there have been great achievements
in improved metrics like accuracy, recall, and F1-score for the de-
identification systems, but also the flexibility of the systems allow the
entities to utilize the same algorithm on different corpora of EMRs. This
property means that an entity can train the system on one corpus of records
and then share that trained system to de-identify another body of reports, or
even share it with other entities. This property, also known as transfer
learning, while very useful and promising, raises the concern of leakage of
sensitive information through these sharing processes (Melis et al. 2019).

One of the most widely used RNN is Long Short-Term Memory (LSTM) architecture,
which is applied in many deep-learning based de-identification systems (Lample
et al. 2016; Kim, Heider, and Meystre 2018; Dernoncourt et al. 2016; Shickel
et al. 2017).
One of these successful algorithms is NeuroNER (Dernoncourt et al. 2016). LSTM
provides a great capability for processing the long dependencies which is an
important characteristic in text formats, thus, providing a very powerful
tool. However, since it has a great number of variables (for instance more
than 40000 for a single LSTM with 100 units), there is the concern of
memorization of data and leakage of the sensitive information in the trained
model due to the model complexity. All these extensive and expensive effort of
the de-identification systems is to protect the privacy of the individuals.
The leak of information through the parameters of the de-identifications model
can jeopardize the privacy and nullify the main goal. One of the recent works
on investigating the data memorization in neural network is (Carlini et al.
2019) on generative sequence models which is a specific neural network with
LSTM. They study these particular LSTMs and show how memorization or leakage
can happen even when the data is rare. Moreover, they suggest a quantitative
metric for measuring the memorization on rare or unique sequences in the data.
While they provide very elegant approach, their method is not directly
applicable to neural networks other than generative ones.

In this article, we focus on the unintended memorization or leakage of the
deep-learning de-identifiers and more specifically, the NeuroNER system. The
reason behind selecting this system is the fact that it has been used in a few
publicly accessible data sets with different structures (i2b2 2014 challenge
data set (Stubbs and Uzuner 2015) and the MIMIC II de-identification data set
(Neamatullah et al. 2008)) and the same system with essentially same structure
and hyperparameters, has been utilized to achieve the state-of-the-art
performance (Dernoncourt, Lee, and Szolovits 2017) (Dernoncourt et al. 2016).
This approach provides a flexible framework and has the potential to be
leveraged in transfer learning paradigms. As such, we are concerned that this
type of approach can encode identities of the individuals in the training data
into the weights of the neural network. In this work we explore this concept.
First we provide a statistical analysis and perform cut-off attacks to
determine the risk if this state-of-the-art algorithm has created unnecessary
exposure for the data subjects. Moreover, we performed membership inference
attack (Shokri et al. 2017), which is the state-of-the-art re-identification
attack related to our approach.

## materials and methods

### NeuroNER

In this subsection, a brief description of the NeuroNER package as well as
specifics of how the system is trained are provided. More details can be found
in (Dernoncourt et al. 2016; Dernoncourt, Lee, and Szolovits 2017). The
statistical and inference attack methods are discussed afterward.

Briefly, the structure of the NeuroNER system, as can be seen in figure 1, is
composed of three layers:

  * •

Token embedding

  * •

LSTM based label prediction

  * •

CRF (conditional random field)

Figure 1: In this chart the circles represent character (Ch) embedding (Emb)
and token level embedding (Embed), the FFNN indicates feed forward neural
network. The attaching points are at the output level of the second layer
after a feed-forward neural network and after the CRF in the third layer as
indicated.

In the first layer, in parallel with a tokenizer and token embedding, there is
an LSTM network for character level embedding. To avoid any memorization or
leak, the token embedding was fixed (no learning) and loaded from the pre-
trained models in public data (as is suggested in the original release codes).
The output of the first layer, which is a concatenation of the word and
character level embedding (illustrated as ’Vec—tor’ in figure 1), is sent to
the bidirectional LSTM of the second layer. The output of the LSTM is then
sent to a feed-forward neural network which will output the probability
vectors $P$. In the third layer, these probability vectors will be the input
of the CRF. The CRF layer can be turned on or off as an option.

The first step in analyzing a text is to break it down to its compartments. In
a sense, one can think of tokenization as separating each word in a text and
calling it token. There are different ready-to-use tokenizers. In this work,
we apply a method known as ‘Spacy’, which is a free and open-source library
for advance natural language processing (NLP).

There are many different methods available to convert tokens into numbers. One
simple approach would be representing any word with one-hot encoding. Given a
vector that has the size of the word domain, the component corresponding to
the word will have value one and all other components in the vector are zero
(so naming is one-hot). In this approach, all the words are perpendicular to
one another and the size of the space grows with the number of words. A more
advanced approach is to use a lower-dimensional (100 for instance) denser
space when the words are represented with vectors that are not all
perpendicular and related words have smaller angles with each other. This
representation is built by an unsupervised learning algorithm that leverages a
large corpus of data, such as Wikipedia, and learns the relevance of the words
in texts (mostly by their co-appearance in the text). In this paper, for the
token embedding, a pre-trained embedding known as Global Vectors for Word
Representation (GloVe) (Pennington, Socher, and Manning 2014) was used which
is of the later type described above.

On the character level, an LSTM with the dimension of 25 is used to embed
tokens into a dense-space using character level training on each token. The
results of the GloVe embedding and character embedding are concatenated as an
input for the next layer.

In NeuroNER, bi-directional LSTM is used in the label prediction in the second
layer.

LSTM is a recurrent neural network architecture and is very effective in
learning from sequence data like text. In the LSTM network in addition to
input and output, each unit has three (input, output and forget) gates
(although there are variations). With the input dimension $m$ and recurrent
dimension $n$, the matrix of weights of input and recurrent connections have a
dimension of $D_{W}=m*n$ and $D_{U}=n*n$, respectively. In addition to the
gates, there is also the state of the unit. There are also biases ($b$) for
each one and so, forward pass one has:

$i^{t}=\sigma(W_{i}x^{t}+U_{i}h^{t-1}+b_{i})$
$f^{t}=\sigma(W_{f}x^{t}+U_{f}h^{t-1}+b_{f})$
$o^{t}=\sigma(W_{o}x^{t}+U_{o}h^{t-1}+b_{o})$
$\tilde{c}^{t}=\tanh(W_{c}x^{t}+U_{c}h^{t-1}+b_{c})$ $c^{t}=f^{t}\circ
c^{t-1}+i^{t}\circ\tilde{c}^{t}$ $h^{t}=o^{t}\circ\tanh(c^{t})$

where $x^{t}(\in\mathbb{R}^{m})$ is the input vector for the unit, $i^{t}$
,$f^{t}$ and $o^{t}$ are input, forget and output activation vectors ,
$\sigma(x)=\frac{1}{1+\exp(-x)}$ and ($\circ$) is element-wise product. The
number of parameters then can be calculated by $4(n*m+n^{2}+n)$. The uni-
directional LSTM, where the sequence can be fed in when the $t$ value
increases. One can make similar calculations when the sequence is trained in
the decreasing order of $t$. In the former case, the network learns from what
comes before a value, while in the later, the network learns from what is
after a value. Bi-directional LSTM is utilizing two LSTMs to learn from all
around the element. This is how the LSTM learns the sequential dependence
between the words.

Conditional random field (CRF) is a powerful statistical method in machine
learning. The conditional part refers to the fact that CRF is a family of
conditional distributions with a structure. For sequence labeling, which is of
interest in this work, linear-chain CRF is the relevant choice which has a
linear structure. In this paper linear-chain CRF is called CRF in short.

### Data sources

The i2b2 challenge data set has been used in this work. The data contains a
set of over 1300 patient records, which is the largest publicly available data
set for de-identification. These reports contain (de-identified) protected
health information (PHI) and it is divided into three sets; namely, training,
validation, and test set. The training set contains more than five hundred
reports.

While the system is trained to label all different PHI types, the patient name
is arguably most sensitive part and there are more than seven hundred patient
surname as well as more than five hundred patient given-name occurrences in
the training data where well over hundreds of them are unique. Here, the focus
is on the patient names.

### Statistical analysis and attack models

To investigate the sharing/leaking of sensitive data, we assume the adversary
has access to almost complete information. Here the almost complete
information simply means that the trained model as well as all the reports are
available to the adversary except that the name which is altered. The goal
here is to investigate the differences in $P$ from data sets that differ in
only one part, names for instance. If the outputs are distinguishable, it
suggests there may be unintentional memorization or leakage of the last names
in the training data. More specifically, new sets of data were produced from
the original one. For the first data set type 1 or inside, names in the
original data set were replaced with other names from the corpus. The next
set, data set type 2 or outside, was constructed by the same procedure but the
original names were replaced with the names from a dictionary of names that
did not appear in the original set. Each of the inside and outside data sets
was divided to three subsets, the one that only surname was replaced (SN), the
one that only given name was replaced (GN) and the one that both surname and
given name were replaced (GN $\&$ SN). The surname dictionary contains above
80000 names, while male and female dictionaries contain about 3000 and 5000
names, respectively. In order to make the comparison with the original
training data, we also preserved punctuation and capitalization of the set and
the names. For statistical attack, the model was trained on original data set
and then used original, data set 1 (inside) and data set 2 (outside) for
testing and compared their relative parameters (the output of softmax, to be
precise).

The non-parametric Kolmogorov—Smirnov (KS) test (Razali, Wah et al. 2011) was
used to determine whether there is any difference in a data set and the
reference probability distribution or between distributions of two data sets
(and thus two-sample KS test). This test gives access to statistic $D$ values.
Statistic $D$ then can be used in numerical or counting algorithms to
estimate/calculate the p-value. To calculate $D$, one has:

$D_{m,n}=\underset{x}{max}|S_{m}(x)-S_{n}(x)|,$

which gives the maximum absolute difference between two distributions with m
and n number of samples. The distributions are relative (normalized) and
empirically produced from samples, and so, sometimes are called empirical
distribution function.

Different re-identification attacks were attempted to stress-check the
vulnerability of the model. They can be categorized to cut-off attacks and
membership inference attack. For the cut-off attacks, the goal was to identify
a set of limits for different probabilities which can be used to re-identify
the participants. In naïve cut-off attack the goal was to find a limit that
can differentiate between the original, data set (1) inside and data set 2
(outside). In brute-force cut-off attack the goal was to feed all possible
names and find if the original name can be re-identified by a limit. For this
purpose, another data preparation was done by randomly selecting three reports
in which the patient name appears six times or more in the body of the report
and then insulating the sentences containing the name (surname) and then
calculating the $P$ for them with replacing the surname with all the over
eighty thousand names in the dictionary.

For membership inference attack, 12 different samples of type data set 2
(outside) were produced (shadow data sets (Shokri et al. 2017)), 10 of which
were used to train 10 different shadow models. As figure 2 illustrates, shadow
model 1 is trained on shadow data 1. Then the training data 1 along with a
randomized mix of reports from shadow data 2 to shadow data 10 (data set -1)
were fed as test to trained model 1 (the mixed data has the same number of
reports as the shadow data 1). The $P$s then with label 1 for shadow data 1
and label -1 for mix shadow data 2 to shadow data 10 were extracted. The same
process was used for all the 10 shadow models. These $P$s and labels then were
fed in to a feed forward network to train it to label inside names (1) and
outside names (-1). The attack accuracy is used as the metric. The shadow
model 11 was used as validation set and shadow model 0 was the target model.
The membership inference attack was conducted where the trained feed-forward
network was used to differentiate and re-identify the names inside shadow data
0 using $P$s produced by testing the shadow model 0 on all possible names
(brute-force) for a few reports with most repetition of a name.

Figure 2: In this chart the ellipse represents data and rectangle indicates
neural network model. The (-k) is the similar data where the names are
replaced with other names that appears in all other sets but set k. Labels are
indicating that the report is in the training data or not, and parameters are
obtained from each model. Figure 3: The graph shows the histogram frequency of
probabilities for the case without CRF. The cyan represents the probabilities
for original/unchanged data. The orange represents the probabilities of
surnames being labeled names when altered randomly with other surnames from
the corpus. The green gives the probabilities for surnames being labeled names
when altered randomly with other surnames from the outside dictionary
(excluding the in-corpus names). The inset zooms in the higher probabilities.
The color-coding is the same as the main graph.

## Parameters and Results

The software was trained on the data set with the following values for the
hyper-parameters: Character embedding dimension and LSTM: 25; Token embedding
dimension and LSTM: 100; Optimizer: SGD; Learning rate: 0.01;Dropout rate:
0.5; Tokenizer: Spacy.

The precision achieved after 95 epoch is 98.43$\%$ on the test set and
99.96$\%$ on the training set. By dropping the third layer (CRF), after 88
epoch the achieved precision is 97.39$\%$ and 98.60$\%$ on the test and
training data respectively. Note that the goal here is to investigate the risk
of leaked information and not to necessarily achieve the best performance on
the test set after training. Nonetheless, the high values of the precision,
recall and F1 indicate that the system has been trained adequately (close to
state-of-the-art) and is reflecting a real-world use of the algorithm. Also,
the suggested precautions, such as freezing the token embedding, were
implemented to maximize the implemented protection methods of the model.

For the inference attack, a feed-forward network with one hidden layer of 64
ReLU units followed by a two softmax units provided the highest accuracy
($\sim$0.75). Please note the inference attack training data is balanced.

The magnitude of the test statistics $D$, and $p$-values obtained from the
double-sided two-sample Kolmogorov-Smirnov tests were used to evaluate the
difference in distributions as shown in table 1 for the network with CRF and
the no-CRF network. As it can be seen, the null hypothesis that the original
names give the same distribution as the altered ones does not hold, especially
in the cases where outside data set 2 (outside) has been used (SN2 and SN2*).
The question remains whether these differences are significant? And more
importantly, can they be utilized for re-identification? Our cut-off attacks
and membership inference attack does not indicate any potential for re-
identification. The naïve cut-off does not indicate any cut-off that can be
used to re-identify any of the reports. The brute-force cut-off attack also
could not narrow down the list of the potential original names to less than
hundreds that would contain the original name used in the training. Even in
the case of membership inference attack, the narrowest of lists got the
original name ranked thirty-eight, while all the rest stayed well above
hundred. In other words, these different attacks failed and no re-
identification was possible.

| no CRF | CRF
---|---|---
| D | p-value | D | p-value
SN1 | 6.2e-2 | 9.5e-2 | 1.0e-1 | 4.3e-4
SN2 | 1.6e-1 | 6.0e-9 | 2.3e-1 | <e-9
SN1* | 6.8e-2 | 5.6e-2 | 9.1e-2 | 3.1e-3
SN2* | 1.7e-1 | <e-9 | 2.3e-1 | <e-9
Table 1: Two-sample Kolmogorov-Smirnov test $D$ statistic and double-sided
p-value comparing the distribution of $P$s for surnames of original
unperturbed data set and other sets which include: SN1 representing the
surnames replaced from the corpus, SN2 surnames replaced from the outside
dictionary, SN1* representing the surnames where both given-names and surnames
have been replaced from the corpus, and SN2* representing surnames where the
surnames and given-names where both replaced from the external dictionaries
(exclusive of names in the corpus).

The histogram (figure 3) illustrates how the distribution of surname
probabilities look like for the random replacement of the original names with
the ones from other reports in the corpus, as well as from outside. The main
graph shows the density of probabilities in different bins for different
surname alterations. The inset figure illustrates the distribution of
probabilities for the highest probabilities only. To the extent that there is
a spread and widening in the statistics, they are overlapping. By using the
cumulative distribution function, some differences can be more easily observed
(figures 4 and 5). Figure 4 shows the cumulative density with Gaussian kernel
density estimates for the case when $P$ is extracted from a model without CRF.
Figure 5 represents the case for $P$ calculated on a network with CRF.

Figure 4: The kernel density estimates are provided for the values of $P$ for
the network with no CRF. Cyan represents the unaltered original data, orange
represents the case when surnames are replaced randomly with other names in
the corpus, green represents the case when the surnames are replaced randomly
with names from outside dictionary(exclusive), red represents $P$ values for
surnames when both given-names and surnames have been replaced from other
names in the corpus, and magenta represents values of $P$ obtained for
surnames when both given-names and surnames were replaced randomly from the
external dictionaries. Figure 5: The kernel density estimates (curves) are
provided for $P$ values for the network with CRF. Cyan represents the
unaltered original data, orange represents the case when surnames are replaced
randomly with other names in the corpus, green represents the case when the
surnames are replaced randomly with names from outside dictionary(exclusive),
red represents $P$ values for surnames when both given-names and surnames have
been replaced from other names in the corpus, and magenta represents $P$
values obtained for surnames when both given-names and surnames were replaced
randomly from the external dictionaries.

## Discussion

In the result section, the difference between distributions for $P$s has been
discussed. These differences are more notable when comparing the original data
set with ones replaced with the external dictionaries. However, one can not
infer that any sensitive information has been compromised. More precisely,
while the p-values are indeed small, one has to note that they have been
gained by sampling over seven hundred surnames. Also, the differences are not
by any means drastic. The overlap of the distributions is overwhelming, and
that can protect sensitive data from adversaries’ inference attacks. For
instance, there is no cutoff that can be used to exclude original $P$s from
either of the other samples.

Moreover, the attempts to narrow down the number of names in the candidates
list, by filtering names with high $P$s between the appearances of the same
patient name in the same report, were also unsuccessful. This failure is due
to the persistence of same names from outside of the corpus with high $P$s,
and also with high fluctuations in the rank of the original name for its
different appearances. Furthermore, it is notable that in this work the almost
absolute knowledge is presumed, meaning that everything including punctuation,
notations, and capitalization of the names were preserved in the replacements
and altered data. Even the membership inference attack did not improve the
case for the leakage of sensitive data as the attack was un-successful to re-
identify any data. Also, the attack could not even limit the candidates for
original name bellow several hundreds as was the case for the cut-off attacks.
To check the robustness of the model, we implemented the membership inference
attack on a reduced data set with fifty reports which is a tenth of the number
of reports in the original set to push the model to over-train and potentially
increase the chance of data leakage, but that did not lead to any leakage.
This indicates all the more to the point that the model does not leak data
when used responsibly. Please note we took suggested precautions as mentioned
in the training section.

It is worth mentioning that $P$ can be interpreted as the probability of
different labels for any input. The distributions of these probabilities for
the original names, as well as replaced names, have very similar properties
like mean, median, standard deviation, and maximum value. That is the case
both for when the model has been trained with CRF layer as well as when CRF
layer has been disabled. In the case of using CRF, the interpretation of the
results of applying softmax on the output of the forward neural network in the
second layer as probabilities is not as accurate. But the numerical values and
statistics of both with and without CRF models, follow each other closely and
the discussions and conclusions provided, stand in both cases.

Differential privacy (Dwork, Roth et al. 2014) (Abadi et al. 2016) has been
gaining momentum in recent years. It is very reliable in the sense that when
applicable, it can provide mathematical insurances to preserve the plausible
deniability up to desired thresholds by introducing noise to the process in a
controlled and measured way. To achieve this mathematically safe-guarded
security, the parameters of this added noise should be set carefully depending
on factors such as the number of epochs, the number of independent data
samples and so on. While new tools like tensor-flow privacy library help with
the implementation of differential privacy, there are caveats for the
practical implementation of the algorithm in cases where the size of available
data is limited. The performance of the algorithms trained in this manner
suffers (Rahman et al. 2018) especially in cases with limited amount of data
and complex models. Moreover, as the number of dependent inputs increases, the
noise should increase, and so the performance gets even worse. Thus, while
theoretically ideal, differential privacy may fall short of enabling
protection of sensitive data by decreasing the accuracy of the model and so
potentially increasing the probability of releasing sensitive data. Of course
it would be interesting to see the above propositions investigated with actual
implementation of differential privacy and its effects, but that is beyond the
scope of this work.

In this paper we have presented an analysis of the potential for a state-of-
the-art deep learning algorithm for de-identification to leak information
concerning subjects’ identities. We show that although statistically there is
difference between the network reaction to inside and outside names, there is
no evidence to suggest that a user could guess that any given subject was
present in the training data, and that the deep neural network encoded the
identities of the users in the training data in a way that creates a risk to
the users. Of course, this does not preclude that some analysis sometime in
the future might reveal the identity of a user, but one can always make the
argument that future technology can do anything, and we feel that this is not
a sufficient argument to be of concern about posting trained networks on
medical data, even when the source training data explicitly contains
identities of individuals. The legal framework governing PHI was developed to
focus on portability (hence the ‘P’ in HIPAA), and judge the trade-off between
risk and benefit of sharing data. This should extend to a new generation of
algorithms trained on such data.

## Acknowledgements

GC, SN and SS are funded by the National Science Foundation, grant # 1822378
‘Leveraging Heterogeneous Data Across International Borders in a Privacy
Preserving Manner for Clinical Deep Learning’. GC and LX are partially
supported by the National Center for Advancing Translational Sciences of the
National Institutes of Health under Award Number UL1TR002378. The content is
solely the responsibility of the authors and does not necessarily represent
the official views of the National Institutes of Health. LX is partially
supported by NIH R01GM118609 Decentralized Differentially-Private Methods for
Dynamic Data Release and Analysis. SN is partially supported by the National
Institutes of Health, award # K01ES025445.

## References

  * Abadi et al. (2016) Abadi, M.; Chu, A.; Goodfellow, I.; McMahan, H. B.; Mironov, I.; Talwar, K.; and Zhang, L. 2016.  Deep learning with differential privacy.  In _Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security_ , 308–318.
  * Carlini et al. (2019) Carlini, N.; Liu, C.; Erlingsson, Ú.; Kos, J.; and Song, D. 2019.  The Secret Sharer: Evaluating and testing unintended memorization in neural networks.  In _28th Security Symposium (Security 19)_ , 267–284.
  * Dernoncourt, Lee, and Szolovits (2017) Dernoncourt, F.; Lee, J. Y.; and Szolovits, P. 2017.  NeuroNER: an easy-to-use program for named-entity recognition based on neural networks.  _Conference on Empirical Methods on Natural Language Processing (EMNLP)_ .
  * Dernoncourt et al. (2016) Dernoncourt, F.; Lee, J. Y.; Uzuner, O.; and Szolovits, P. 2016.  De-identification of Patient Notes with Recurrent Neural Networks.  _Journal of the American Medical Informatics Association (JAMIA)_ .
  * Dwork, Roth et al. (2014) Dwork, C.; Roth, A.; et al. 2014.  The algorithmic foundations of differential privacy.  _Foundations and Trends® in Theoretical Computer Science_ 9(3–4): 211–407.
  * Kim, Heider, and Meystre (2018) Kim, Y.; Heider, P.; and Meystre, S. 2018.  Ensemble-based Methods to Improve De-identification of Electronic Health Record Narratives.  In _AMIA Annual Symposium Proceedings_ , volume 2018, 663. American Medical Informatics Association.
  * Lample et al. (2016) Lample, G.; Ballesteros, M.; Subramanian, S.; Kawakami, K.; and Dyer, C. 2016.  Neural architectures for named entity recognition.  _arXiv preprint arXiv:1603.01360_ .
  * Melis et al. (2019) Melis, L.; Song, C.; Cristofaro, E. D.; and Shmatikov, V. 2019.  Exploiting Unintended Feature Leakage in Collaborative Learning.  In _2019 2019 IEEE Symposium on Security and Privacy (SP)_ , 497–512. Los Alamitos, CA, USA: IEEE Computer Society.  ISSN 2375-1207.  doi:10.1109/SP.2019.00029.
  * Neamatullah et al. (2008) Neamatullah, I.; Douglass, M. M.; Lehman, L.-w. H.; Reisner, A.; Villarroel, M.; Long, W. J.; Szolovits, P.; Moody, G. B.; Mark, R. G.; and Clifford, G. D. 2008.  Automated de-identification of free-text medical records.  _BMC Medical Informatics and Decision Making_ 8(1): 32.  ISSN 1472-6947.  doi:10.1186/1472-6947-8-32.
  * Pennington, Socher, and Manning (2014) Pennington, J.; Socher, R.; and Manning, C. D. 2014.  GloVe: Global Vectors for Word Representation.  In _Empirical Methods in Natural Language Processing (EMNLP)_ , 1532–1543.
  * Rahman et al. (2018) Rahman, M. A.; Rahman, T.; Laganière, R.; Mohammed, N.; and Wang, Y. 2018.  Membership Inference Attack against Differentially Private Deep Learning Model.  _Trans. Data Priv._ 11(1): 61–79.
  * Razali, Wah et al. (2011) Razali, N. M.; Wah, Y. B.; et al. 2011.  Power comparisons of shapiro-wilk, kolmogorov-smirnov, lilliefors and anderson-darling tests.  _Journal of Statistical Modeling and Analytics_ 2(1): 21–33.
  * Shickel et al. (2017) Shickel, B.; Tighe, P. J.; Bihorac, A.; and Rashidi, P. 2017.  Deep EHR: a survey of recent advances in deep learning techniques for electronic health record (EHR) analysis.  _IEEE Journal of Biomedical and Health Informatics_ 22(5): 1589–1604.
  * Shokri et al. (2017) Shokri, R.; Stronati, M.; Song, C.; and Shmatikov, V. 2017.  Membership inference attacks against machine learning models.  In _2017 IEEE Symposium on Security and Privacy (SP)_ , 3–18. IEEE.
  * Stubbs and Uzuner (2015) Stubbs, A.; and Uzuner, Ö. 2015.  Annotating longitudinal clinical narratives for de-identification: The 2014 i2b2/UTHealth corpus.  _Journal of Biomedical Informatics_ 58: S20–S29.
  * Yang and Garibaldi (2015) Yang, H.; and Garibaldi, J. M. 2015.  Automatic detection of protected health information from clinic narratives.  _Journal of Biomedical Informatics_ 58: S30–S38.


# LiDAR Based Detection and Classification of Pedestrians and Vehicles Using
Machine Learning Methods

Farzad Shafiei Dizaji
Engineering Systems and Environment
University of Virginia
Charlottesville, VA 22904
<EMAIL_ADDRESS>

###### Abstract

The goal of this paper is to classify objects mapped by LiDAR sensor into
different classes such as vehicles, pedestrians and bikers. Utilizing a LiDAR-
based object detector and Neural Networks-based classifier, a novel real-time
object detection is presented essentially with respect to aid self-driving
vehicles in recognizing and classifying other objects encountered in the
course of driving and proceed accordingly. We discuss our work using machine
learning methods to tackle a common high-level problem found in machine
learning applications for self-driving cars: the classification of pointcloud
data obtained from a 3D LiDAR sensor.

_K_ eywords LiDAR data, Machine Learning, Self-driving

##  1 Introduction

LiDAR is a surveying approach that measures distance to a target by
illuminating that target with a pulsed laser light and measuring the reflected
pulses with a sensor. Differences in laser return times and wavelength can
then be used to make digital 3D representations of the target [1].

There has been considerable work on object classification for autonomous
vehicle navigation. Most of the works use camera vision to capture 2D images
of vehicles and other vulnerable road users such as pedestrians and
bicyclists. However, camera vision alone may not be able to provide important
depth information to detect and track objects with the level of reliability
needed for safe driving. LiDAR enables generation of 3D images from a single
sweep. Multiple sweeps can be used to generate information about velocities
and distances. For greater accuracy, we may use multi-sensor fused data, which
is already used for traffic light detection [2]. In our project, we intend to
train the data using Neural Networks after any required reprocessing or
feature extraction. Neural Networks use the processing of the brain as a basis
to develop algorithms that can be used to model complex patterns or prediction
problems [3]. The basic network architecture of neural networks has an input
layer, hidden layer(s) and an output layer. The hidden layers are used as
‘filters’ that make the network faster and efficient by identifying the
important patterns of input that must be passed on to the next layers. Neural
networks have the ability to learn from initial inputs and their
relationships, generalize the model and predict on unseen data. Our goal is to
create a system that can efficiently and accurately predict and classify ahead
of time the objects in its surroundings, so as to aid the self-driving car in
maneuvering through safely.

The LiDAR sensor data is present as a dense three-dimensional point cloud.
Before it can be used for classification, a lot of pre-processing needs to be
done to refine the data, cluster the point cloud and reduce the number of
clusters. To begin with, the point cloud frame data is extracted, following
which ground filtering is carried out in an effort to facilitate the
clustering process. The Random Sample Consensus (RANSAC) algorithm is used for
the ground filtering process. RANSAC operates by producing best fit linear or
quadratic planes by observing objects in the environments as outliers. From
there, we can isolate grounded points by thresholding points in the z-axis.
Next, “Density Based Spatial Clustering of Application with Noise” (DBSCAN)
algorithm is used for clustering. It is a density-based algorithm which is
simple and fast and does not require any initialization and allows for an
unlimited number of clusters. The approach produces clusters based on
euclidean distance criteria and a minimum number of samples. Following
clustering, valid and usable features are extracted from the clusters. This
generates useful candidates that can be fed as input data to machine learning
algorithms for classification via Support Vector Machines and Artificial
Neural Networks.

As autonomous vehicles advance towards handling realistic road traffic, they
face street scenarios where the dynamics of other traffic participants must be
considered explicitly [4]. LiDAR enables extensive 3D mapping of its
environment, which can be used to navigate a self-driving car or robot
predictably and safely. By using LiDAR sensor data, it is possible to generate
position data for objects surrounding the car ahead of time [1, 5]. This
becomes essential in achieving full autonomy in self driving cars, especially
in areas with unpredictable traffic situations. In this project, we seek to
enable quick classification of LiDAR sensor data into vehicles, pedestrians,
bikers, etc. to enable a self-driving car to plan its trajectory accordingly.

The paper is presented as follows: we formally introduce our problem in
section 2 and provide some background in section 3. We then discuss the KITTI
dataset and the preprocessing needed to perform feature extraction and
classification operations in sections 4, 5. We then describe our
classification architecture and our experimental results in 6. Finally, we
conclude with information on our software and hardware and some final remarks
in section 7 and section 8.

##  2 Problem Definition

Consider an autonomous agent in an unknown environment with no extra
information about the potential shape, dynamics, or behaviors of the
environment. Our agent seeks to identify objects in this environment given by
structural data captured in the form of pointclouds which are defined as
tuples containing an x, y, z coordinate and an intensity value. Given a
pointcloud O, there exists some subset of O that includes points belonging to
a specified class. In our work, we consider three classes: pedestrians, cars,
and cyclists. We desire then to correctly label each point so that it belongs
to a subset Oped, Ocar, Ocyclist, and Oignored (to indicate points which fit
no class but can be safely ignored such as walls or trees).

##  3 Background

In existing literature, we can identify a plethora of work with a large surge
since the 2007 DARPA Urban challenge. Dr. Sebastian Thrun at Stanford
University has been a pivotal investigator of this work, producing a great
wealth of research on object tracking using shape, color, and motion [6],
object tracking with semisupervised learning [7], and segmentation of RGBD
data [8]. From these contributions, a standard method for handling 3D
pointcloud data has emerged. First, if there are multiple sources of LiDAR
data, the pointclouds are fused according to some registration technique. In
many cases, this is conducted via Iterative Closest Point (ICP), an
optimization algorithm which seeks to minimize the error between corresponding
points in two pointclouds. From there, the data is segmented into objects of
interest. Segmentation draws a parallel in machine learning terminology in
clustering as its major goal is to produce groupings of interesting data
points which may be used for tracking and classification. Segmentation is an
active research topic but is commonly performed using one of a variety of
techniques including DBSCAN, neural network classification, and geometric
fitting. Being real-time is one of the most important characteristics of an
algorithm in this field. But there are numerous factors that can impact this
capability such as path planning. A rapid path planning algorithm is needed
that can generate the shortest path from an initial point to a goal point to
gather the essential information from the surrounding environment. Recently an
exact geometry based path planning algorithm has been presented in [9, 10]
which has $O(nn^{\prime 2})$ running time where $n$ shows the number of the
vertices and $n^{\prime}$ is the number of effective polygons which can be
convex or non-convex. This algorithm has shown its capabilities in finding a
collision-free and shortest path among a group of polygons in a very short
time. Also, in order to produce real-time speeds, the technical space becomes
more limited and algorithms like DBSCAN are used less due to their long
processing times. Finally, after the objects have been segmented, the data can
be evaluated in parallel for tracking and classification. Tracking extracts
relevant state information on the cluster subject such as velocity, volume,
position, reachable space, etc. This is commonly performed using an Extended
Kalman Filter that attempts to figure out the corresponding A, B, C, and D
matrices in state-space convention that dictate the dynamics of the object.
Classification attempts to take relevant information about the object and
identify what type of object it is, such as a cyclist, pedestrian or a car.
Generally, these do not need the information extracted from tracking but new
research is utilizing some of the motion characteristics of an object to help
classification. This is commonly performed with some feature extraction and
then either artificial neural networks or support vector machines. In this
work, we choose to show our results in testing a SVM, ANN, and a decision
tree.

Figure 1: Sample label in KITTI dataset with two cars and a cyclist

##  4 Dataset

We use the KITTI dataset [11], to train our algorithms. The KITTI dataset,
released in 2012 contains several hours of LiDAR footage with object labels
for various classes including our desired pedestrian, cyclist, and car
classes. The dataset is also widely recognized as an excellent method for
standardized testing as the dataset appears in many works of literature. In
Figure 1, a sample set of labeled points is shown clustered by color. In the
image, a cyclist and two vehicles can be somewhat clearly identified.

###  4.1 Data and Preprocessing

We currently plan to use the KITTI Vision dataset, which is publicly available
and contains six hours worth of frame-by-frame traffic footage. The data
contains both camera images and LiDAR point maps, as well as tags for which
images contain which objects [12]. We hope to use these features and
compatibility with existing software to which we have access to make it
simpler to classify the data.

Depending on the format and labeling of the data, we may need to perform our
own registration and segmentation on the data. If given multiple streams of
LiDAR data, we may need to fuse these different sources from their different
perspectives to achieve better results. This can be performed easily using the
open source PCL library and its variety of registration algorithms such as
iterative closest point or feature matching through the use of SIFT features
[13]. Segmentation takes the raw data and finds continuities, grouping
densities of data into classifiable objects. The segmented data can then be
used for machine learning applications [14].

###  4.2 Algorithms

The algorithms that we intend to use for classification in our project are
Neural Networks and Support Vector Machines. Neural Networks are greatly used
in pattern recognition and classification problems, owing to their ability to
take in a lot of inputs, and process them to infer complex, non-linear
relationships [3]. Support Vector Machines are classifiers that use labeled
training data to generate an optimal hyper-plane that can be used to classify
unseen data [15, 16].

##  5 Pre-processing

Since the KITTI datset presents completely raw data, we must utilize a way to
filter out ground data and to cluster the objects appropriately. These are
both sources of error in our work but also exist as new potential avenues of
future research and have appropriately large amounts of literature. In our
case, we chose to use a method of ground filtering produced by an author at
Beijing Normal University known as Cloth Simulation Filtering (CSF) [17] and
the Mean Shift algorithm for clustering.

###  5.1 CSF Ground Filtering

CSF ground filtration relies on a simulated model of a cloth commonly utilized
by computer graphics. The cloth is assigned a certain rigidity and is laid on
an inverted version of the point-cloud. Due to the natural behaviors of the
cloth, rolling hills in the cloth are captured while sharp changes indicating
objects are ignored. The method requires some tuning but yielded excellent
results. As a side note, assigning extreme rigidity to the cloth essentially
forces this algorithm to become a plane fitting algorithm.

###  5.2 Mean-Shift Clustering

Mean-shift clustering is an extension of Kernel Density Estimation (KDE), a
technique used to estimate an underlying probability density function of a
specific set of data points. Mean-shift clustering uses a Gaussian kernel in
KDE to evaluate a data-set and observe strong groupings in the data. Each
point then is iteratively moved to the closest grouping of data, creating
natural clusters based on densities of points. A nice advantage of this
approach is that the approach does not require knowledge of the number of
clusters to identify which is an important for our work. The approach however,
also suffers from speed in larger frames of point-cloud data although for our
approach the speed was sufficient (Figure 2).

Figure 2: Ground Filtering via CSF - Ground Points are marked in Yellow &
Clustering via Mean-Shift - Clusters are marked by different colors

##  6 Pipeline

Features were selected based on features found in existing literature. In many
works we found common reference to the use of eigenvalues or eigenvectors to
indicate the directionality of a pointcloud. For example, a pedestrian would
have a large eigenvalue in the upward direction since the characteristic shape
of a pedestrian is thinner than a cyclist or a vehicle. In the end, we
selected five features found in literature which were the following. The
eigenvalues of the covariance matrix for each cluster in x, y, and z
directions, the total volume of the object, and the variance in intensity of
the cluster. These features were extracted from 23 of the testing data
clusters identified by the object labels for training while the remaining 13
were used for validation.

Because the KITTI dataset was very large, we found that it was extremely
difficult to utilize k-fold cross validation for our approach due to the
limited amount of data we had available and the limited amount of
computational time we had to conduct an adequately size k-fold cross
validation. Since a single extraction and training sequence took several hours
due to the ground filtering and clustering pre-processing, we decided that
breaking the dataset into folds was possible but entirely too time consuming
given our current hardware. However, we discuss here that model selection
approach is easily doable with enough time and strong enough computing power.
In total, we were able to extract features from 16,314 car samples, 1,016
cyclist samples, and 2,918 pedestrian samples.

Our experimental design utilized the previous section’s work to train and
validate a support vector machine and a decision tree. The approach is
graphically shown in Figure 3. After conducting the pre-processing stage of
the pointcloud data, we fit/train our algorithms using the features extracted
from the clusters of each frame of the training set. Then, the data is used to
classify the validation data for testing purposes to identify quality.

Figure 3: Experimental Architecture
Figure 4: Experimental Architecture

##  7 Results

We present the results of our work with the following metrics. Due to the
unique nature of our dataset, we present the data based on the percentages of
each frame that was correctly annotated. In general, we found that our
approach yielded a highly noisy but accurate classification methodology. The
metrics are presented as the total percentage of data points correctly
classified and the total percentage of the labeled data points correctly
classified. The reasoning behind these differing metrics is due to the extreme
density of points per frame. Since a significant portion of points belong to
no such class, such as walls or trees, the first metric shows how the
algorithm should theoretically only classify points in the frame that are of
interest. The second metric shows that of the points that are actually of
interest, how well does the algorithm perform. For the decision tree and
support vector machines, we yielded an overwhelming number of false positives
but excellent accuracy. In other words, we found that both approaches yielded
high recall but low precision. Our results are presented in the following
(Figure 4).

There were a few large sources of error that were easy to understand after
going through testing. The first was that we identified objects as ignorable
when our classifier had a low prediction score, for our case we used below 90
percent. However, without specifying a class for do not cares, we found that
most time our algorithm would just attempt to classify everything with a high
margin. This resulted in extremely noisy data and produced incredibly
difficult results to understand. Another source of error were for pedestrians
and pole shaped objects. Tree trunks, signs, and telephone poles were often
misclassified as pedestrians due to their similar eigenvalues and volumes. In
general, it would be safe to say that even the intensity variances were
comparable. Finally, a final source of error was that our clustering algorithm
struggled with objects that were beyond a certain size (Figure 5). When tuning
the clustering parameters, it became difficult to select a bandwidth that was
thin enough that close objects were correctly segmented while wide enough that
larger objects were not segmented into several components. In future attempts,
we could introduce a new layer of pre-processing that may attempt to fuse
multiple nearby clusters if there were enough continuous points (Table 1).

Multi-layer neural networks which are sometimes referred to as perceptrons,
are simple models of several connected neurons similar to what is seen in the
natural neural networks of animals. The main objective of building these
artificial models of brain has been to design systems which can show some
learning capabilities like the natural brain [18, 19]. We also attempted to
classify LiDAR data using neural networks. We used MLFFNNs to build our model
and will use cross-validation to check the fit of the model (Table 2). After
studying the dataset, it is clear that the data is skewed and will require
proper resampling to train the model effectively. Our plan is to resample the
data such that we have a balanced dataset, train our model and check the fit
of our model for solving the problem and to choose the best parameters to use
in our model. Using the 5 previously mentioned features, we chose only to
classify pedestrians and vehicles with deep perceptrons.

Table 1: Number of states based on team parameters Accuracy
---
Classifier | Total Frame | Total Cluster
SVM | 0.04 | 0.638
DT | 0.04 | 0.676

Figure 5: Misclassification Examples

According to the data, Class 1 (pedestrian) is under-represented and accounts
for only 3 % of the whole dataset. If we train our model using this dataset,
the model will be inefficient and it will be trained to predict only Class 0
because it will not have sufficient training data. One of the other problems
of using this skewed dataset to train our model is that since Class 1 is
under-represented, the model will assume that it is a rare case and will try
to predict positive due to the lack of training data. We may get a high
accuracy when we test our model but we should not be confused by this because
our dataset does not have a balanced test data. Hence, we have to rely on the
recall which relies on TP and FP. In cases where we have skewed data, adding
additional data of the under-represented feature (over-sampling) is an option.
Since we don’t have that option we have to resort to under-sampling. Under-
sampling of the dataset involves keeping all our under-represented data (Class
1) while adding the same number of features of Class 0 (cars) to create a new
dataset comprising of an equal representation from both classes. So, after
applying Under-sampling method, we have to train/test our model using 2032
rows of data. It is not much but it will have to suffice. Our approach was
conducted via MATLAB, sci-kit Learn, and Keras using a standard i7 Intel CPU.
Future implementations and approaches could be heavily improved through the
use of a GPU and parallel code as the Mean-Shift algorithm is extremely easy
to run in parallel.

##  8 Conclusion

great deal of pre-processing to filter ground data and extract important
features for our classifiers. These elements greatly impacted our approach
speed, presenting nowhere near real-time speeds. This is reflected in existing
literature as new literature is beginning to move away from SVM and feature-
extraction methods and move towards comprehensive neural network based
approaches which simultaneously perform feature extraction and classification
extremely quickly. These approaches are also naturally useful since they can
be leveraged for all aspects of the processing pipeline from classification to
segmentation for clustering and for ground filtering. However, as neural
networks become more prevalent in self-driving literature an important
question arises in how these methodologies work. Many researchers criticize
ANN’s for their black-box nature, being able to produce amazing results with
no knowledge of their internal workings. This has an extreme relevance in our
future culture as future issues that arise with neural networks may be far
more difficult to resolve and may lead to catastrophic failures in autonomous
systems that rely on their capabilities.

Table 2: Table to test captions and labels Layer(type) | Activation Function | Output Shape | Parameters
---|---|---|---
Dense1 | Dense3 | (None, 200) | 1200
Dropout1 | n/a | (None, 200) | 0
Dense2 | relu | (None, 200) | 40200
Dropout2 | n/a | (None, 200) | 0
Dense3 | Sigmoid | (None, 1) | 201
Table 3: Loss Architectures | Training Acc | Training Loss | Validation Acc | Validation Loss
---|---|---|---|---
Dense Layer 2 | 0.93 | 0.16 | 0.92 | 0.17
Dense Layer 3 | 0.95 | 0.16 | 0.95 | 0.17
Dense Layer 4 | 0.95 | 0.13 | 0.95 | 0.14

## References

  * [1] C. Premebida, G. Monteiro, U. Nunes, and P. Peixoto, “A lidar and vision-based approach for pedestrian and vehicle detection and tracking,” in _2007 IEEE intelligent transportation systems conference_. IEEE, 2007, pp. 1044–1049.
  * [2] N. Fairfield and C. Urmson, “Traffic light mapping and detection,” in _2011 IEEE International Conference on Robotics and Automation_. IEEE, 2011, pp. 5421–5426.
  * [3] S. Lawrence, C. L. Giles, A. C. Tsoi, and A. D. Back, “Face recognition: A convolutional neural-network approach,” _IEEE transactions on neural networks_ , vol. 8, no. 1, pp. 98–113, 1997.
  * [4] J. Levinson, J. Askeland, J. Becker, J. Dolson, D. Held, S. Kammel, J. Z. Kolter, D. Langer, O. Pink, V. Pratt _et al._ , “Towards fully autonomous driving: Systems and algorithms,” in _2011 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2011, pp. 163–168.
  * [5] K. Kidono, T. Miyasaka, A. Watanabe, T. Naito, and J. Miura, “Pedestrian recognition using high-definition lidar,” in _2011 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2011, pp. 405–410.
  * [6] D. Held, J. Levinson, S. Thrun, and S. Savarese, “Combining 3d shape, color, and motion for robust anytime tracking.” in _Robotics: science and systems_ , 2014.
  * [7] A. Teichman and S. Thrun, “Tracking-based semi-supervised learning,” _The International Journal of Robotics Research_ , vol. 31, no. 7, pp. 804–818, 2012.
  * [8] A. Teichman, J. T. Lussier, and S. Thrun, “Learning to segment and track in rgbd,” _IEEE Transactions on Automation Science and Engineering_ , vol. 10, no. 4, pp. 841–852, 2013.
  * [9] H. Jafarzadeh, S. Gholami, and R. Bashirzadeh, “A new effective algorithm for on-line robot motion planning,” _Decision Science Letters_ , vol. 3, no. 1, pp. 121–130, 2014.
  * [10] H. Jafarzadeh and C. H. Fleming, “An exact geometry–based algorithm for path planning,” _International Journal of Applied Mathematics and Computer Science_ , vol. 28, no. 3, pp. 493–504, 2018.
  * [11] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for autonomous driving? the kitti vision benchmark suite,” in _2012 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 2012, pp. 3354–3361.
  * [12] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The kitti dataset,” _The International Journal of Robotics Research_ , vol. 32, no. 11, pp. 1231–1237, 2013.
  * [13] R. B. Rusu and S. Cousins, “Point cloud library (pcl),” in _2011 IEEE international conference on robotics and automation_ , 2011, pp. 1–4.
  * [14] E. Grilli, F. Menna, and F. Remondino, “A review of point clouds segmentation and classification algorithms,” _The International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences_ , vol. 42, p. 339, 2017.
  * [15] Y. Sang, H. Zhang, and L. Zuo, “Least squares support vector machine classifiers using pcnns,” in _2008 IEEE Conference on Cybernetics and Intelligent Systems_. IEEE, 2008, pp. 290–295.
  * [16] N. Cristianini, J. Shawe-Taylor _et al._ , _An introduction to support vector machines and other kernel-based learning methods_. Cambridge university press, 2000.
  * [17] W. Zhang, J. Qi, P. Wan, H. Wang, D. Xie, X. Wang, and G. Yan, “An easy-to-use airborne lidar data filtering method based on cloth simulation,” _Remote Sensing_ , vol. 8, no. 6, p. 501, 2016.
  * [18] A. Joghataie and M. S. Dizaji, “Nonlinear analysis of concrete gravity dams by neural networks,” in _Proceedings of the World Congress on Engineering_ , vol. 2. Citeseer, 2009.
  * [19] A. Joghataie, M. S. Dizaji, and F. S. Dizaji, “Neural network software for dam-reservoir foundation interaction,” in _Proceedings of the International Conference on Mechanical, Automotive and Materials Engineering, Dubai, United Arab Emirates_ , 2012, pp. 7–8.


# ICO learning as a measure of transient chaos in $\mathcal{PT}$-symmetric
Liénard systems

J. P. Deka<EMAIL_ADDRESS>Department of Physics,
Girijananda University, Guwahati - 781017, India  A. Govindarajan
<EMAIL_ADDRESS>Department of Nonlinear Dynamics, School of Physics,
Bharathidasan University, Tiruchirappalli - 620024, India  A. K. Sarma
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology
Guwahati, Guwahati - 781039, India

###### Abstract

In this article, we investigate the implications of the unsupervised learning
rule known as Input-Correlations (ICO) learning in the nonlinear dynamics of
two linearly coupled $\mathcal{PT}$-symmetric Liénard oscillators. The fixed
points of the oscillator have been evaluated analytically and the Jacobian
linearization is employed to study their stability. We find that on increasing
the amplitude of the external periodic drive, the system exhibits period-
doubling cascade to chaos within a specific parametric regime wherein we
observe emergent chaotic dynamics. We further notice that the system indicates
an intermittency route to chaos in the chaotic regime. Finally, in the
period-4 regime of our bifurcation analysis, we predict the emergence of
transient chaos which eventually settles down to a period-2 oscillator
response which has been further validated by both the maximal Finite-Time
Lyapunov Exponent ($\mathcal{FTLE}$) using the well-known Gram-Schmidt
orthogonalization technique and the Hilbert Transform of the time-series. In
the transiently chaotic regime, we deploy the ICO learning to analyze the
time-series from which we identify that when the chaotic evolution transforms
into periodic dynamics, the synaptic weight associated with the time-series of
the loss oscillator exhibits stationary temporal evolution. This signifies
that in the periodic regime, there is no overlap between the filtered signals
obtained from the time-series of the coupled $\mathcal{PT}$-symmetric
oscillators. In addition, the temporal evolution of the weight associated with
the stimulus mimics the behaviour of the Hilbert transform of the time-series.

##  I Introduction

In neuroscience, Hebbian learning is an attempt to explain the synaptic
efficacy when the post-synaptic cells are continuously stimulated by the pre-
synaptic cells [1]. In this context, it could be quoted that ‘Neurons that
fire together wire together’ [4]. Synaptic plasticity plays an important role
in the coordinated implementation of cells that fire in synchrony.
Furthermore, there is an evidence that synaptic strength could either be
strengthened or weakened due to heterosynaptic plasticity. In this context,
the ICO learning rule is a method that incorporates the rules of classical
Hebbian learning (CHL) along with the conditions of heterosynaptic plasticity
[1-3,5-6, 44]. It takes into account the correlation between the synapses and
based on this, the synaptic strength (measured in terms of the weight) could
be strengthened or weakened. In addition, this rule is a variant of
differential Hebbian learning and it takes into account the changes in the
post-synaptic cells. One major advantage of this method is that it circumvents
the potentially destabilising term in CHL by replacing the derivative of the
output of the learning rule with the derivative of the input. This enables the
use of higher learning rates together with the removal of the destabilising
term in the learning rule.

On the other hand, Carl M. Bender and his graduate student Stefan Boettcher’s
ground-breaking discovery [7] that non-Hermitian Hamiltonians can yield a real
eigenspectra paved the path for a new era in the foundational studies of
quantum mechanics. Such Hamiltonians were termed as $\mathcal{PT}$-symmetric
Hamiltonians and as such, they are invariant under the joint operation of
parity ($\mathcal{P}$) and time-reversal ($\mathcal{T}$) operators by
commuting with the $\mathcal{PT}$-operator $[H,PT]=0$. Unlike Hermitian
Hamiltonians, such non-Hermitian Hamiltonians possess an
exceptional-point ($\mathcal{EP}$) which is accompanied by the so-called
$\mathcal{PT}$-symmetric phase transition in the eigenspectra after which the
former changes from real to imaginary. The regime where the eigenspectrum is
real is known as unbroken $\mathcal{PT}$-symmetric regime and if the
eigenspectrum is purely imaginary, then the regime is called as the broken
$\mathcal{PT}$-symmetry.

In the context of experimental endeavors, El-Ganainy et al. [8] were the first
to put forward the proposition that evanescently coupled dielectric optical
slab-based waveguides with balanced gain and loss can facilitate a feasible
way to validate these mathematical theories in optics. Such a proposition was
deemed conceivable due to the fact that the time-independent Schrödinger
equation and the paraxial equation of diffraction/dispersion are isomorphic.
The potential function in the former and the refractive index in the latter
serve as the connecting tool for this shared isomorphic trait. Not long ago,
the first experimental observation was reported in a coupled optical waveguide
configuration by Rüter et al. [8] followed by Guo et al. [9]. Two optical
waveguides were etched on to a lithium niobate substrate. One of the channels
was optically pumped to provide amplification while the other waveguide was
provided an equal proportion of attenuation. In such a
$\mathcal{PT}$-symmetric optical waveguide configuration, the $\mathcal{PT}$
phase transition was characterized by a transition from periodically evolving
optical power to exponential growth and decay of optical power in the two
channels. In the context of optical systems, the notion of
$\mathcal{PT}$-symmetry has been explored in diverse domains such as
optoelectronic oscillators [10], optomechanical systems [11-13], modulational
instability in complex media [14], active LRC circuits [15], multilayered
structures [16-17], unidirectional invisibility in periodic structures
[18-19], wireless power transfer [20], optical lattices [21-24], double ring
resonators [25], and so on.

The temporal dynamics of most continuous nonlinear systems are usually either
asymptotically stable (otherwise oscillatory) or unstable. Oscillatory
temporal dynamics, on the other hand, could either be periodic, quasiperiodic
or chaotic. But in transient chaotic dynamics, the temporal dynamics undergo a
transition in their behavior from chaotic to periodic as the system is allowed
to evolve. Hence, it can be inferred that chaotic dynamics that arise for a
finite length of time can be termed as transient chaos [26-31]. Grebogi, Ott
and Yorke’s discovery of the emergence of transient chaotic dynamics in
certain systems passing through crisis led to the birth of this field in
nonlinear dynamics [29]. Moreover, Cruthfield and Kaneko’s work on the
interconnection between attractors and turbulence (spatiotemporal complexity)
generated a new perspective on the study of turbulence. Such phenomena are
related to attractors or long-lived transient dynamics of non-attracting
chaotic sets [30]. Subsequently, transient chaos has been reported in
electronic oscillators [32-34], lasers [35-37], parametrically forced pendulum
[38], Duffing oscillators [39] and so on. Quantification of such chaotic
dynamics requires the computation of the maximal $\mathcal{FTLE}$, which has
been observed to decrease from a positive value to zero, thereby signifying
the emergence of a limit cycle attractor [39].

In this article, we report the transient chaotic dynamics that emerge in a
$\mathcal{PT}$-symmetric coupled Liénard systems for non-autonomous
excitation. The article is organized as follows. In section II, we present a
brief overview of the mathematical model. This is followed by an analysis of
the stability of the fixed points using the Jacobian linearization and how
transient chaotic dynamics arise in our system, which we have further
validated by studying the maximal $\mathcal{FTLE}$ and the Hilbert transform
of the time-series in section III. And finally, in Section IV, we investigate
the implication of the unsupervised learning rule known as ICO learning in the
analysis of time-series transformation from chaotic to peridic evolution,
followed by our conclusion in Sec. V.

##  II MATHEMATICAL MODEL

Let us consider the $\mathcal{PT}$-symmetric Liénard oscillator given below.

$\ddot{x}+f(x)\dot{x}+\alpha x+g(x)=0$ (1)

where $\alpha$ is a real constant, $f(x)$ is the dissipation/amplification
profile and $g(x)$ is the nonlinearity profile of the oscillator. Under the
joint operation of the $\mathcal{PT}$-operator ($x\rightarrow-x$ and
$t\rightarrow-t$), one can show that the above equation gets transformed as
follows

$\ddot{x}-f(-x)\dot{x}+\alpha x-g(-x)=0$ (2)

Now, if $f(x)$ and $g(x)$ are both anti-symmetric in $x$, then the oscillator
whose dynamics is governed by Eq. (1) is $\mathcal{PT}$-symmetric. So, in this
article, we consider two configurations of coupled $\mathcal{PT}$-symmetric
oscillators with dissipation profile $f(x_{i})=\eta x_{i}$ and nonlinearity
profile $g(x_{i})=\beta x_{i}^{3}$ (where $i=1,2$) given below

$\displaystyle\ddot{x}_{1}+\eta x_{1}\dot{x}_{1}+\alpha x_{1}+\beta
x_{1}^{3}+\kappa x_{2}=0$ (3a) $\displaystyle\ddot{x}_{2}+\eta
x_{2}\dot{x}_{2}+\alpha x_{2}+\beta x_{2}^{3}+\kappa x_{1}=0$ (3b)

Here, $\kappa$ is the coupling constant. Also, the term $\eta$ denotes the
strength of nonlinear dissipation/amplification and $\beta$ refers to the
coefficient of Duffing nonlinearity of the two oscillators. It may be noted
here that the dissipation/amplification profile $f(x_{i})=\eta x_{i}$ could
provide both dissipation and amplification in the same oscillator depending on
whether the parameter $\eta$ is positive or negative. It must also be noted
that in this autonomous case, the individual oscillators are $\mathcal{PT}$
symmetric, but the coupled oscillator configuration, as a whole, is not
$\mathcal{PT}$ symmetric.

In the context of an external agent driving the system, the mathematical model
will be transformed as follows.

$\displaystyle\ddot{x}_{1}-\eta x_{1}\dot{x}_{1}+\alpha x_{1}+\beta
x_{1}^{3}+\kappa x_{2}=A_{0}cos(\omega t)$ (4a)
$\displaystyle\ddot{x}_{2}+\eta x_{2}\dot{x}_{2}+\alpha x_{2}+\beta
x_{2}^{3}+\kappa x_{1}=A_{0}cos(\omega t)$ (4b)

where $A_{0}$ is the amplitude of the external drive and $\omega$ is the
frequency of the drive. One can easily observe that this system remains
invariant under the joint operation of the parity operator defined as
$\mathcal{P}:x_{1}\leftrightarrow x_{2}$ and the time-reversal operator
defined as $\mathcal{T}:t\rightarrow-t$ and hence, it is obvious that this
coupled oscillator configuration is $\mathcal{PT}$-symmetric. It must also be
noted that the choice of the external drive is vital in such a scenario and as
such, we have chosen a cosine function drive which is an even function of
time. Now, since the amplification/dissipation profile chosen in both
configurations is similar, in the second configuration, the term signifying
gain/loss in one of the oscillators should carry the opposite sign so as to
play the antithetical role. Hence, we will term the first system as the gain
and the second as the loss oscillator in our results and discussion section.
In our analysis, to excite the gain oscillator we have chosen the initial
condition as, $(x_{1},y_{1},x_{2},y_{2})=(0.1,0,0,0)$, where
$y_{i}=\dot{x}_{i}$ and $i=1,2$.

##  III ANALYSIS OF FIXED POINTS AND TIME-SERIES

Before analyzing the coupled oscillator configuration, we would like to study
the stability of the fixed points of the first oscillator in the absence of
coupling for the autonomous case. By setting $\dot{x}_{1}=y_{1}$ and
$\kappa=0$, we have

$\displaystyle\dot{x}_{1}=y_{1}$ (5a) $\displaystyle\dot{y}_{1}=-\eta
x_{1}y_{1}-\alpha x_{1}-\beta x_{1}^{3}$ (5b)

Hence the fixed points (FP) of this nonlinear system are read as

  1. 1.

FP1 $\rightarrow$ $(x_{1},y_{1})=(0,0)$

  2. 2.

FP2 $\rightarrow$ $(x_{1},y_{1})=(\sqrt{-\alpha/\beta},0)$

  3. 3.

FP3 $\rightarrow$ $(x_{1},y_{1})=(-\sqrt{-\alpha/\beta},0)$,

Figure 1: Eigenspectra of the Jacobian linearization of (a) FP1 (b) FP2 and
(c) FP3 for $\beta=0.5$ and $\alpha=-1$

Spectral analysis of the Jacobian eigenspectra shows that FP1 is a saddle
fixed point, FP2 is a stable fixed point and FP3 is an unstable fixed point,
which have been depicted clearly in Fig. 1. One interesting feature that one
can observe from our spectral analysis is that our system contains all three
forms of hyperbolic fixed points. Moreover, it is to be noted that the
potential function for a single oscillator is given by $V(x)=\alpha
x_{1}^{2}/2+\beta x_{1}^{4}/4$ and stability analysis of the equilibrium
points of this potential attests that the equilibrium point
$x=\pm\sqrt{-\alpha/\beta}$ is stable if $\alpha<0$ and $\beta>0$. Generally,
we have a double well potential for certain choice of parameters. But the
presence of nonlinear position dependent dissipation renders one of the stable
equilibrium points unstable. Furthermore, another interesting aspect to be
noted here is that the eigenspectra of all three fixed points exhibit no
exceptional point (EP). Normally, in a $\mathcal{PT}$-symmetric system, the
eigenspectra undergo a phase transition from real to imaginary eigenvalues as
certain parameters (like gain/loss) are varied. But in this configuration
where the position dependent dissipation ($\eta$) makes the system to be
$\mathcal{PT}$-symmetric, there is no such a phase transition in the
eigenspectra of the Jacobian linearization albeit the system obeys the so-
called $\mathcal{PT}$-symmetry. A complete and detailed analysis of the
eigenspectra for the linearization Jacobian for this system has been presented
in our previous work [43].

We would now like to start our discussion with the bifurcation diagram of the
temporal maxima of the system under study versus the amplitude of the external
drive $A_{0}$. It can be clearly seen in Fig. 2 that on increasing the value
of $A_{0}$, the temporal dynamics of the gain oscillator transforms from
periodic to chaotic. We have considered a very lengthy time-series for our
bifurcation analysis, ranging from $t=0$ to $t=25000$, out of which we have
removed the temporal data up to $t=23000$ as transient. We will now
concentrate on a few specific cases of the amplitude of the external drive.
Some of the regimes in the gain oscillator response to the amplitude of the
external drive are shown in Table I. The phase plane of the oscillator
response for this parametric regime is shown in Fig. 3. In our analysis, we
have chosen the following parameter values: $\eta=0.1$, $\beta=0.5$,
$\kappa=0.5$, $\alpha=0.25$ and $\omega=0.25$. The amplitude of the external
drive $A_{0}$ is chosen as the variable parameter.

Figure 2: Temporal maxima of the gain-oscillator $x_{1,max}$ vs. $A_{0}$.
Table 1: Temporal response of the gain oscillator.

It can be clearly seen from Fig. 4 that as we increase the value of $A_{0}$,
the gain oscillator starts to exhibit multiple temporal maxima and beyond a
certain threshold, we observe the emergence of chaotic dynamics in the
oscillators. This has also been portrayed through the phase plane analysis in
Fig. 3. We would now like to present an in-depth analysis of the parametric
regime before the system transits into chaotic dynamics. One can notice from
Fig. 4 that the system exhibits a case of period doubling cascade to chaos.
But at the same time, the bifurcation analysis also depicts chaotic dynamics
for $A_{0}=0.27148$ and $A_{0}=0.27772$ and a period-2 and a period-6
oscillator response for $A_{0}=0.2784$ and $A_{0}=0.27852$, respectively. This
means even within the periodic regime, there are regions where the oscillators
might portray chaotic or period-N temporal dynamics ($N\geq 2$). This has been
depicted in the time-series of the gain oscillator shown in Fig. 5. Hence, we
can draw a conclusion that our system exhibits intermittency in the
bifurcation analysis depicting the period doubling cascade to chaos.

Figure 3: Phase plane of the gain-oscillator for (a) $A_{0}=0.17$, (b)
$A_{0}=0.205$, (c) $A_{0}=0.212$ and (d) $A_{0}=0.246$. Figure 4: Temporal
maxima of the gain-oscillator delineating the period doubling cascade to
chaos. Figure 5: Temporal evolution of the gain-oscillator for (a)
$A_{0}=0.27148$, (b) $A_{0}=0.27772$, (c) $A_{0}=0.2784$ and (d)
$A_{0}=0.27852$.

Also, from Fig. 6, we can observe that the temporal evolution delineate
intermittent chaotic dynamics for $A_{0}=0.282$. This could be attributed to
the reason that $x_{1}=-0.7$ is a stable fixed point for the gain oscillator
and this has been shown in the spectral analysis of the Jacobian linearization
in our previous work [43]. So, from the temporal evolution in Fig. 6, one can
notice that for the duration between $t=1240$ to $t=1440$, the gain oscillator
exhibits a period-4 oscillatory dynamics about its stable fixed point. Hence,
we can conclude that $\mathcal{PT}$-symmetric Liénard systems reveal the
intermittency route to chaos. Furthermore, intermittency could be easily
characterized by studying the temporal maxima of the time-series. In fact, it
is the same algorithm that we used for our bifurcation analysis in Fig. 2 and
Fig. 4. This has been plotted in Fig. 7 from which one can notice that there
is a tiny window from $t=1240$ to $t=1440$ wherein we observe that there is a
periodicity in the temporal evolution of $x_{1,max}$.

Figure 6: Temporal evolution of the gain-oscillator for $A_{0}=0.282$. Figure
7: Temporal maxima of the gain-oscillator for $A_{0}=0.282$.

Now, we would like to concentrate on a particular case in the parametric
regime of Fig. 4. From Fig. 8(a), one can easily identify that for
$A_{0}=0.274$, the system exhibits chaotic dynamics before it settles down
into periodic sinusoidal oscillations. One of the peculiar characteristics
that can be perceived from Fig. 8(a) is that the system appears to be in a
transiently chaotic state for an incredibly long time. This, in fact, proved
quite a tedious computational challenge when the temporal maxima of the time-
series were being evaluated for the bifurcation diagrams in Fig. 2 and Fig. 4.
This is the reason why we have considered an extremely long time-series for
our bifurcation analysis. From Fig. 8(b) and 8(c), it could be gathered that
the chaotic phase plane eventually transforms into a period-2 temporal
evolution. Transient chaos can be easily characterized by analyzing the
maximal $\mathcal{FTLE}$. As the temporal evolution of the two oscillators is
allowed to progress, it is observed that the chaotic oscillations eventually
cease after a finite length of time and that the time-series cease to be
chaotic after $t\approx 5000$ (Fig. 8(a)). This leads us to infer that our
system exhibits transient chaotic dynamics. In other words, it means that
chaos exists for a finite length of time. Transient chaos can be characterized
by studying the maximal $\mathcal{FTLE}$ with a moving window of finite time-
series. We have chosen the Gram-Schmidt orthogonalization method to calculate
the maximal $\mathcal{FTLE}$ [40-42]. So, accordingly, in our numerical
computation, we have chosen to evaluate the $\mathcal{FTLE}$ for a
sufficiently long temporal evolution of the system. Considering these aspects,
we have plotted Fig. 8(d) which depicts the maximal $\mathcal{FTLE}$ against
the length of the time-series. As observed in the figure, while we increase
the length of $t$, $\mathcal{FTLE}$ exponentially decays and there is an
asymptotic convergence of $\mathcal{FTLE}$ towards the $x$-axis. At the same
time, if we look at the time-series in Fig. 8(a), one can notice that the
chaotic oscillations cease after $t=5000$. Upon further increasing the value
of $t$, the asymptotic convergence of $\mathcal{FTLE}$ continues which
validates that the oscillations try to converge to a limit cycle attractor.
Hence, from our analysis, one can claim that when we have a coupled
$\mathcal{PT}$-symmetric Liénard system with the same characteristics except
for the nonlinear dissipation/amplification and the system is driven by an
external agent, the oscillators exhibit the emergence of a chaotic attractor
which ceases to exist after a finite length of time as the oscillating
temporal dynamics converge to a limit cycle attractor (Fig. 8(c)).

Figure 8: (a) Time-series, (b) transient state phase plane (c) steady-state
phase-plane and (d) $\mathcal{FTLE}$ of the gain oscillator for $A_{0}=0.274$.

Transient chaos can further be characterized by studying the Hilbert Transform
of the time-series which is given by the following relation,

$\hat{H}[x_{i}(t)]=\frac{1}{\pi}P.V.\int_{-\infty}^{+\infty}\frac{x_{i}(\tau)}{t-\tau}d\tau$
(6)

where $x_{i}(t)$ is the time-series of the system. In Fig. 9, we have plotted
the Hilbert Transform of the time-series shown in Fig. 8(a). From the figure,
one can observe that in the regime where the time-series is transiently
chaotic, the Hilbert Transform is seen to increase monotonically up to
$t\approx 5000$ after which the time-series transforms into a periodic
temporal evolution. At this junction, the Hilbert Transform yields a constant
value as the oscillator exhibits periodic temporal evolution. Thus, we confirm
that studying the Hilbert Transform of the time-series serves as a viable way
to characterize transient chaos in the nonlinear systems [43].

Figure 9: Hilbert transform of the time-series of the gain oscillator for
$A_{0}=0.274$.

##  IV ICO Learning

We shall now explore the implications of the unsupervised learning rule known
as ICO Learning in our transiently chaotic time-series of the
$\mathcal{PT}$-symmetric coupled oscillator system. This rule was first put
forward by Prof. Florentin Wörgötter and Prof. Bernd Porr, which is a
modification of Hebbian differential learning [1-3, 44]. In 2023, the optics
group under Prof. Claudio R. Mirasso and Prof. Ingo Fischer of IFISC-UIB,
Spain experimentally demonstrated this technique using a pulsed-optoelectronic
oscillator [45]. In this rule, the synaptic weights account for the
instantaneous change in presynaptic connection by the cross-correlation
between the presynaptic events instead of the presynaptic and postsynaptic
events. Here, there are two types of input signals - the single reference
signal $u_{1}$ and the multiple stimulus signals $u_{i>1}$ and their
corresponding weights are updated by the rule given below.

$\displaystyle\frac{dw_{i}}{dt}=\eta u_{i}\frac{du_{1}}{dt}$ (7a)

where $\eta=0.01$ is the learning rate in our simulations. The weight
associated with the reference signal is fixed at $w_{1}=1$ and the weight for
the stimulus $w_{i>1}$ evolves with time and $u_{1}$ and $u_{i>1}$ are the
temporally extended spiking signals. This rule is an instance of unsupervised
neural algorithm and it measures the temporal separation between successive
neural events. If we have one reference and one stimulus signal, then when
$u_{1}$ precedes the $u_{2}$, the weight $w_{2}$ decreases and when the
reverse happens, $w_{2}$ increases. This is a signature of heterosynaptic
plasticity, where the evolution of synaptic strength are triggered by the
activity in other pathways.

Figure 10: (a-b) Temporal maxima and (c-d) filtered temporal maxima of the
time-series of the gain and loss oscillators for $A_{0}=0.274$. Figure 11:
Temporal evolution of the weight $w_{2}$ in the temporal regime (a)
$3.961\times 10^{5}$ to $3.965\times 10^{5}$ and (b) $3.965\times 10^{5}$ to
$3.968\times 10^{5}$. Other parameters are, $Q=0.51$, $c=1$ and $f=0.01$.
Figure 12: (a) Temporal evolution of the weights $w_{1}$ and $w_{2}$ and (b)
the output $v$ for $A_{0}=0.274$.

In our system, we have two oscillators with continuous-time evolution. We
consider the time-series of the gain oscillator $x_{1}$ as the reference
signal and that of the loss oscillator $x_{2}$ as the stimulus signal. Now,
the first task is to evaluate the temporal maxima in the time-series, where we
used the technique employed in plotting the bifurcation diagram in Fig. 2 and
Fig. 4. In addition, we shall consider the time-series of transiently chaotic
dynamics of our system in discrete time-steps. If a temporal maximum was
evaluated, we consider that as $1$ and if no temporal maximum was evaluated,
we take that as $0$ in the time-series for temporal maxima. Afterwards, to
increase the temporal overlap in the spiking signals, the raw spiking signal
$x_{i,max}$ (which contains the data of the temporal maxima) is filtered and
the function given below is employed.

$\displaystyle h(n)=\frac{1}{dc}e^{an}sin(dn)$ (8a)

where $n$ is the discrete time-step, $a=-\pi f/Q$ and $d=\sqrt{(2\pi
f)^{2}-a^{2}}$. Here, $c>0$ is a parameter to tune the amplitude of the
filtered signal, $f$ is frequency of the filtered normalized to unit sampling
rate and $Q>0.5$ is the decay rate. Our configuration is a two oscillator
system and so, we shall use two filters for the raw signals containing
information of the spiking dynamics to transform them into the activation
signals $u_{1}$ and $u_{2}$. The output of this learning procedure will be
generated as a linear superposition of the filtered signals as follows

$\displaystyle v=w_{1}u_{1}+w_{2}u_{2}$ (9a)

Fig. 10 depicts two scenarios in the evolution of the temporal maxima of the
two signals. In Fig. 10(a), we can see that the stimulus (red) precedes the
reference (black) and in Fig. 10(b), the reference precedes the stimulus. The
same could also be seen for the filtered signals in Figs. 10(c) and 10(d).
Under both scenarios, it could be seen from Fig. 11(a) that when the stimulus
$u_{2}$ precedes the reference $u_{1}$ as shown in Fig. 10(c), $w_{2}$ could
be seen to increasing in a spiking way and then, it decreases. But overall,
$w_{2}$ could be seen to be increasing. But when $u_{1}$ precedes $u_{2}$ as
shown in Fig. 10(d), $w_{2}$ is seen to be decreasing in Fig. 11(b).
Furthermore, such spiking is behavior is observed at the time-step when the
filtered signals overlap each other. From these facts, we can thus ascertain
that our simulations provide us with the correct result. Now in Fig. 12(a),
when we run our simulation till the $8\times 10^{5}$ time-step, the temporal
evolution of $w_{2}$ is erratic but on an average, it increases in the chaotic
regime. But when the system starts exhibit periodic dynamics, $w_{2}$ becomes
stationary. This means that the filtered signals of the reference and stimulus
no longer overlap and the temporal separation between the two is constant.
Likewise, the output $v$ also becomes periodic with no intermittent spiking as
observed in the chaotic regime in Fig. 12(b). Note that a similar result was
also observed in Fig. 9 where we have evaluated the Hilbert Transform of the
time-series of the gain oscillator. All these findings clearly suggest that
both the ICO Learning rule and Hilbert Transform serve a similar approach in
the analysis of transiently chaotic dynamics in nonlinear systems.

##  V Conclusion

In conclusion, we have witnessed the period-doubling cascade to chaos with a
brief window in which we observed intermittent chaotic behaviour in the
system. The oscillators in this model also display transient chaotic dynamics
in a certain parametric regimes, as validated by the maximal $\mathcal{FTLE}$.
In addition, we observed that the Hilbert transform of the time-series exhibit
a transformation from erratic to stationary as the time-series transform from
chaotic to periodic. Also, we have applied the unsupervised learning rule
called the ICO learning to the analysis of our time-series. When the system
exhibits chaotic dynamics before transforming into periodic, there is an
intermittent increment as well as decrement in the weight associated with the
stimulus which means that the temporal separation between the signals changes
erratically with time. But when the time-series transforms from chaotic to
periodic, the weight associated with the stimulus becomes stationary
signifying no overlap between the reference and the stimulus filtered signals.
Hence, from our analysis, one can argue that this technique is an efficient
measure in the analysis of the temporal evolution of any system and its
associated transformation from aperiodic to periodic both in the experimental
and computational sectors. Furthermore, it could be said that this learning
procedure plays a role similar to the Hilbert transform in time-series
analysis.

## ACKNOWLEDGEMENTS

J.P.D would like to acknowledge his group members - Prof. Claudio R. Mirasso,
Prof. Ingo Fischer, Dr. Miguel C. Soriano, Dr. Apostolos Argyris and Dr.
Silvia Ortín Gonzalez at the Instituto de Fisica Interdisciplinar y Sistemas
Complejos (IFISC) of the Universitat de les Illes Balears (UIB), Spain for the
valuable research experience in the context of ICO Learning.

  1. [1]

D. O. Hebb, The organization of behavior; a neuropsychological theory, Wiley
(1949).

  2. [2]

B. Porr, F. Wörgötter, Strongly Improved Stability and Faster Convergence of
Temporal Sequence Learning by Using Input Correlations Only, Neural
Computation 18 (6), 1380–1412 (2006).

  3. [3]

K. Möller, D. Kappel, M. Tamosiunaite, C. Tetzlaff, B. Porr, F. Wörgötter,
Differential Hebbian learning with time-continuous signals for active noise
reduction, PLoS ONE 17, 5 (2022).

  4. [4]

S. Löwel, W. Singer, Selection of intrinsic horizontal connections in the
visual cortex by correlated neuronal activity, Science 255, 209–12 (1992).

  5. [5]

S. Royer, D. Paré, Conservation of total synaptic weight through balanced
synaptic depression and potentiation, Nature 422, 518 (2003).

  6. [6]

C. H. Bailey, M. Giustetto, Y. Y. Huang, R. D. Hawkins, E. R. Kandel, Is
heterosynaptic modulation essential for stabilizing hebbian plasiticity and
memory, Nat. Rev. Neurosci. 1, 11 (2000).

  7. [7]

C. M. Bender and S. Boettcher, Real Spectra in Non-Hermitian Hamiltonians
Having $\mathcal{PT}$ Symmetry, Phys. Rev. Lett. 80, 5243 (1998).

  8. [8]

C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, D.
Kip, Observation of parity-time symmetry in optics, Nat. Phys. 6, 192-195
(2010).

  9. [9]

A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez,
G. A. Siviloglou, and D. N. Christodoulides, Observation of
$\mathcal{PT}$-Symmetry Breaking in Complex Optical Potentials, Phys. Rev.
Lett. 103, 093902 (2009).

  10. [10]

J. Zhang and J. Yao, Parity-time–symmetric optoelectronic oscillator, Science
Advances 4, 6 (2018).

  11. [11]

X. Lü, Hui Jing, Jin-Yong Ma and Ying Wu, $\mathcal{PT}$-Symmetry-Breaking
Chaos in Optomechanics, Phys. Rev. Lett. 114, 253601 (2015).

  12. [12]

W. Li, Y. Jiang, C. Li and H. Song, Parity-time-symmetry enhanced
optomechanically-induced-transparency, Sci. Rep. 6, 31095 (2016).

  13. [13]

X. Xu, Y. Liu, C. Sun, and Y. Li, Mechanical $\mathcal{PT}$ symmetry in
coupled optomechanical systems, Phys. Rev. A 92, 013852 (2015).

  14. [14]

A. K. Sarma, Modulation instability in nonlinear complex parity-time symmetric
periodic structures, J. Opt. Soc. Am. B 31, 18611866 (2014).

  15. [15]

J. Schindler, A. Li, M. C. Zheng, F. M. Ellis, and T. Kottos, Experimental
study of active LRC circuits with $\mathcal{PT}$ symmetries, Phys. Rev. A 84,
040101(R) (2011).

  16. [16]

M. Sarisaman, Unidirectional reflectionlessness and invisibility in the TE and
TM modes of a $\mathcal{PT}$-symmetric slab system, Phys. Rev. A 95, 013806
(2017).

  17. [17]

J. P. Deka and A. K. Sarma, Highly amplified light transmission in a parity-
time symmetric multilayered structure, Appl. Opt. 57, 1119 (2018).

  18. [18]

Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, and D. N.
Christodoulides, Unidirectional Invisibility Induced by
$\mathcal{PT}$-Symmetric Periodic Structures, Phys. Rev. Lett. 106, 213901
(2011).

  19. [19]

L. Feng, Y. Xu, W. S. Fegadolli, M. Lu, J. E. B. Oliveira, V. R. Almeida, Y.
Chen and A. Scherer, Experimental demonstration of a unidirectional
reflectionless parity-time metamaterial at optical frequencies, Nat. Mater.
12, 108-113 (2013).

  20. [20]

S. Assawaworrarit, X. Yu and S. Fan, Robust and efficient wireless power
transfer using a switch-mode implementation of a nonlinear parity–time
symmetric circuit, Nat. 546, 387 (2017).

  21. [21]

A. Regensburger, M. -Ali Miri, C. Bersch, J. Näger, G. Onishchukov, D. N.
Christodoulides, and U. Peschel, Observation of Defect States in PT-Symmetric
Optical Lattices, Phys. Rev. Lett. 110, 223902 (2013).

  22. [22]

A. Regensburger, C. Bersch, M. -Ali Miri, G. Onishchukov, D. N.
Christodoulides and U. Peschel, Parity–time synthetic photonic lattices, Nat.
488, 167 (2012).

  23. [23]

M. -Ali Miri, A. Regensburger, U. Peschel, and D. N. Christodoulides, Optical
mesh lattices with $\mathcal{PT}$ symmetry, Phys. Rev. A 86, 023807 (2012).

  24. [24]

M. Wimmer, A. Regensburger, M. -Ali Miri, C. Bersch, D. N. Christodoulides and
U. Peschel, Observation of optical solitons in $\mathcal{PT}$-symmetric
lattices, Nat. Comm. 6, 7782 (2015).

  25. [25]

J. P. Deka and A. K. Sarma, Chaotic dynamics and optical power saturation in
parity–time ($\mathcal{PT}$) symmetric double-ring resonator, Nonlinear Dyn.
96, 565 (2019).

  26. [26]

T. Tél, The joy of transient chaos, Chaos 25, 097619 (2015).

  27. [27]

Y. -C. Lai and T. Tél, Transient Chaos: Complex Dynamics on Finite-Time Scales
(Springer, New York, 2011).

  28. [28]

C. Grebogi, E. Ott, and J. Yorke, Chaotic Attractors in Crisis, Phys. Rev.
Lett. 48, 1507 (1982).

  29. [29]

C. Grebogi, E. Ott, and J. Yorke, Crises, sudden changes in chaotic
attractors, and transient chaos, Physica D 7, 181 (1983).

  30. [30]

C. Dettmann and G. Morriss, Hamiltonian reformulation and pairing of Lyapunov
exponents for Nosé-Hoover dynamics, Phys. Rev. E 55, 3693 (1997).

  31. [31]

J. R. Crutchfield and K. Kaneko, Are Attractors Relevant to Turbulence?, Phys.
Rev. Lett. 60, 2715 (1988).

  32. [32]

F. T. Arecchi and F. Lisi, Hopping Mechanism Generating 1/f Noise in Nonlinear
Systems, Phys. Rev. Lett. 49, 94 (1982).

  33. [33]

R. Hilborn, Quantitative measurement of the parameter dependence of the onset
of a crisis in a driven nonlinear oscillator, Phys. Rev. A 31, 378 (1985).

  34. [34]

R. Rollins and E. Hunt, Intermittent transient chaos at interior crises in the
diode resonator, Phys. Rev. A 29, 3327 (1984).

  35. [35]

F. T. Arecchi, R. Meucci, G. Puccioni, and J. Tredicce, Experimental Evidence
of Subharmonic Bifurcations, Multistability, and Turbulence in a Q-Switched
Gas Laser, Phys. Rev. Lett. 49, 1217 (1982).

  36. [36]

D. Dangoisse, P. Glorieux, and D. Hennequin, Laser Chaotic Attractors in
Crisis, Phys. Rev. Lett. 57, 2657 (1986).

  37. [37]

F. Papoff, D. Dangoisse, E. Poite-Hanoteau, and P. Glorieux, Chaotic
transients in a $CO_{2}$ laser with modulated parameters: Critical slowing-
down and crisis-induced intermittency, Opt. Commun. 67, 358 (1988).

  38. [38]

R. Leven, B. Pompe, C. Wilke, and B. Koch, Determining Lyapunov exponents from
a time series, Physica D 16, 285 (1985).

  39. [39]

S. Sabarathinam, K. Thamilmaran, Transient chaos in a globally coupled system
of nearly conservative Hamiltonian Duffing oscillators, Chaos, Solitons and
Fractals 73, 129 (2015).

  40. [40]

A. Uchida, Optical Communication with Chaotic Lasers: Applications of
Nonlinear Dynamics and Synchronization (John Wiley $\&$ Sons, 2012).

  41. [41]

F. Christiansen and H. H. Rugh, Computing Lyapunov spectra with continuous
Gram - Schmidt orthonormalization, Nonlinearity 10, 1063 (1997).

  42. [42]

I. S. Gomez, Lyapunov exponents and poles in a non Hermitian dynamics, Chaos,
Solitons and Fractals 99, 155 (2017).

  43. [43]

J. P. Deka, A. K. Sarma, A. Govindarajan, M. Kulkarni, Multifaceted nonlinear
dynamics in $\mathcal{PT}$-symmetric coupled Liénard oscillators, Nonlinear
Dyn 100, 1629 (2020).

  44. [44]

B. Porr, F. Wörgötter, Isotropic sequence order learning, Neural Comput 15,
831 (2003).

  45. [45]

S. Ortín, M. C. Soriano, C. Tetzlaff, F. Wörgötter F, I. Fischer, C. R.
Mirasso and A. Argyris, Implementation of input correlation learning with an
optoelectronic dendritic unit, Front. Phys. 11, 1112295 (2023).