$\ast$$\ast$affiliationtext: Equal contribution.11affiliationtext: Warsaw
University of Technology

# Rethinking of Encoder-based Warm-start Methods in Hyperparameter
Optimization

Dawid Płudowski Antoni Zajko Anna Kozak Katarzyna Woźnica

###### Abstract

Effectively representing heterogeneous tabular datasets for meta-learning
remains an open problem. Previous approaches rely on predefined meta-features,
for example, statistical measures or landmarkers. Encoder-based models, such
as Dataset2Vec, allow us to extract significant meta-features automatically
without human intervention. This research introduces a novel encoder-based
representation of tabular datasets implemented within the liltab package
available on GitHub 111https://github.com/azoz01/liltab. Our package is based
on an established model for heterogeneous tabular data proposed in [Iwata and
Kumagai, 2020]. The proposed approach employs a different model for encoding
feature relationships, generating alternative representations compared to
existing methods like Dataset2Vec. Both of them leverage the fundamental
assumption of dataset similarity learning. In this work, we evaluate
Dataset2Vec and liltab on two common meta-tasks – representing entire datasets
and hyperparameter optimization warm-start. However, validation on an
independent metaMIMIC dataset highlights the nuanced challenges in
representation learning. We show that general representations may not suffice
for some meta-tasks where requirements are not explicitly considered during
extraction.

##  1 Introduction

Due to their inherent heterogeneity, tabular datasets pose a unique challenge
in the realm of data representation [Vanschoren, 2019, Jomaa et al., 2021].
Unlike image or text data, they exhibit diverse variables, dimensions, and the
absence of intrinsic variable order. Despite these complexities, tabular data
remains ubiquitous across various domains [Davenport and Kalakota, 2019,
Alanazi, 2022]. From that perspective, this is relevant to develop their
effective representations suitable for meta-learning and knowledge transfer
methods [Vanschoren, 2019].

So far, existing approaches to encoding diverse tabular datasets into common
vector space primarily relied on predefined meta-features. Some aggregate
information about the statistics of one-dimensional variables contained in the
datasets and some are based on information theory or landmarkers –
characteristics of specific machine learning models[Rivolli et al., 2022].
Still, these definitions were proposed in advance and mostly based on human
intuition and available information extraction methods. Yet, at the core of
their definitions lay a crucial oversight: they do not explicitly consider how
the resulting representations would be applied within the context of meta-
learning applications.

In the case of tabular datasets, meta-learning is narrowed to extracting
information on how to build models for new, unseen data. The unsolved problem
is primarily the problem of algorithm selection and HP optimization of
classical machine learning models. Most often, the quality of the
representation of tabular data sets was evaluated in terms of accelerating
Bayesian Optimization (BO) [Snoek et al., 2012] by using warm-start determined
from the most similar datasets. Traditional a priori defined meta methods have
shown some improvement in terms of speed up of the optimization [Wistuba et
al., 2015a, b, 2018, Feurer et al., 2018]. However, it is possible that better
representations can be found that have a higher correlation with the task for
which they are intended to be employed.

The emergence of encoders tailored to heterogeneous datasets opens new
possibilities for representation learning. The main advantage over
representations composed of predefined meta-features is the ability to
formalize desired properties of the representations within the loss function.
Pioneering encoder-based representations of heterogeneous tabular datasets in
Dataset2Vec extractor, Jomaa et al. [2021] introduce:

######  Requirement 1.

Ensuring close representations for batches of observations within the same
dataset while maintaining distinct representations for observations across
different datasets.

By imposing this condition, the authors assume that the resulting
representations will reflect the internal complex structure of the data. Then,
the representations are used to classify whether a set of observations is
derived from the same dataset. In this particular task, the representations
prove to be effective. The meta-task of warm-starting BO was verified as an
additional one. However, the broader impact of the encoder-based
representations remains open. Furthermore, the lack of another meta-feature
extractor hinders the further verification of the potential of encoder-based
representations in a defined meta-task.

To fill these gaps, we propose a novel approach to extract representations
from the encoder component of the model introduced in [Iwata and Kumagai,
2020]. Selected architecture is designed explicitly for heterogeneous
datasets. We employ the fundamental assumption of dataset similarity learning
and introduce new dataset representations implemented within the liltab
package. Our approach, while rooted in Requirement 1, benefits from a
different encoder-based model to capture feature relationships, thereby
generating representations alternative to Dataset2Vec.

To evaluate the efficacy of liltab representations, we conduct a comprehensive
analysis using purposeful designed datasets from metaMIMIC [Woźnica et al.,
2023] along with well-established dataset repositories such as OpenML
[Vanschoren et al., 2013] and UCI [Markelle Kelly, 2021]. Our primary meta-
task focuses on optimizing hyperparameters of classical machine learning
models, serving as a litmus test for representation quality. Precisely, we
assess the efficiency of warm-start points in hyperparameter optimization
through Bayesian Optimization.

The validation process highlights the nuanced challenges in representation
learning. Merely adhering to Requirement 1 through encoder-based
representations proves inadequate for broader meta-learning.

##  2 liltab’s representations

We implemented the Python package – liltab, which provides an interface to the
predictive model proposed in [Iwata and Kumagai, 2020] and extends its
functionalities to an encoding of the entire dataset. The unique feature of
Iwata’s model is the ability to perform predictions on any tabular data in the
supervised setting. The predictions are obtained using an inference network
that casts batches of heterogeneous datasets to a common latent space.

Liltab uses latent vectors for features and targets representations from the
re-implemented model. We focus on the sub-component of the predictive
architecture – the inference network – where latent representations of
specific observations are obtained. Achieving such encoding is done in three
steps. In the first step, the model extracts marginal distributions of the
features and targets by propagating each value in the dataset through a neural
network, followed by pooling. In the second step, it concatenates obtained
vectors, passes them through the second neural network, and aggregates its
outputs. This operation results in a fixed-size vector for each observation,
which implicitly catches joint distributions. Finally, it uses the last neural
network to cast them into the final latent space of observations’
representations. To ensure that dataset representation does not depend on the
number of observations, we average all vectors corresponding to specific
records.

During the original training procedure, the authors of [Iwata and Kumagai,
2020] wanted to maximize its predictive performance. However, representations
produced by an inference network occur to be uninformative on real-world data
and often are vectors with few non-zero coordinates. To overcome that issue,
we propose a new encoder-based technique that targets this part of the model.

In our training procedure, each batch of data is created from samples from
several different datasets, where the numbers of samples for each dataset are
equal. Next, observations from one dataset in the batch are encoded to their
representation. During the calculation of the loss function, we focus on
representations of specific observations and train the network to produce
encodings that are close to each other when they correspond to the same
dataset and far away when they relate to different datasets. To achieve that,
we propose the contrastive learning approach with loss function inspired by
[Song et al., 2015]:

$loss=\frac{1}{|\mathcal{N}|+|\mathcal{S}|}\left(\sum_{i,j\in\mathcal{S}}\|f(x_{i})-f(x_{j})\|_{2}-\sum_{i,j\in\mathcal{N}}\|f(x_{i})-f(x_{j})\|_{2}\right)\text{,}$

where $x_{i}$ denotes specific observation, $f$ is a liltab’s encoder
function, $\mathcal{S}$ is a set of pairs of indices of observations that
belong to the same dataset, and $\mathcal{N}$ denotes pairs of indices of
observations that are from different datasets. Note that we index by
observations, not by subsets in batch.

##  3 Considered approaches

Next to the liltab-based encoder, we use Dataset2Vec encoder meeting
Requirement 1. Both work on heterogeneous data, which is crucial when working
with tabular data. Their performance is contrasted with two baseline methods
independent of dataset representations.

###  3.1 Dataset2Vec’s representations

Dataset2Vec (D2V) is an encoder of the whole dataset proposed in [Jomaa et
al., 2021]. It is one of the approaches to encoding heterogeneous tabular data
that uses neural networks instead of statistical measures. It is a DeepSet-
based method [Zaheer et al., 2017] composed of feed-forward neural networks.
The main idea of this model is to treat datasets as a hierarchical set, i.e.,
a set of features and targets contains sets of values. Firstly, the model
captures an interdependency between features and corresponding targets. Next,
it pools output from the previous layer and aggregates across all features and
targets. Finally, the model projects the output into meta-features space,
resulting in dataset encoding, a single vector with constant dimensionality.
Authors of D2V show that these representations satisfy the Requirement 1,
which is crucial for our evaluation.

###  3.2 Baselines

As baselines for the usability of encoder-based representations, we use two
methods: BO without warm-start and warm-start based on the ranks of HP
configurations. We define ranks by the number of tasks in which the HP
configuration was the best. In the case of metaMIMIC tasks, because of the
small size of the tasks, many HP configurations occur to result in the same
best score for specific tasks, which introduces noise to the ranks. To prevent
such issues, we consider the configuration the best one when the number of the
best HP configurations is at most three in a specific HP optimization task.
Otherwise, no HP configuration is taken into account in determining the ranks.

###  3.3 Methods comparison

Here, we want to emphasize the differences between the two presented encoders:
D2V and liltab’s inference network.

Although they both produce vectors as the representation of the dataset, they
have significantly different architectures; specifically, D2V was designed to
be a datasets classifier while a predictive model, that is a core for liltab
encoder, is used to perform predictions on heterogeneous data in few-shot
learning manner.

The preferred size of the data also differs. D2V is trained, both in the
original paper and in our experiments, on medium-size data (i.e., hundreds or
thousands of records), and liltab uses its few-shot learning architecture.

Finally, there is a subtle distinction between optimization functions in the
training of both encoders: for liltab, the optimization ends on the level of
single observations while D2V is optimized on the level of the whole subsets
of the dataset.

Regarding the baselines, we want to emphasize the need to use both BO without
warm-start and a method based on ranks. The first approach is the simplest
method that does not use a database with scores of HP configurations. The
second one is also a simple method, yet uses the mentioned database. In some
cases, one or the other baseline may perform better as the HP database can be
both an advantage of the algorithm and its limitation. When the knowledge
collected on different tasks is not transferable to the new task, methods
based on the database may be misled and start optimization in points which
lead to suboptimal solutions. In such cases, a baseline without a warm-start
will perform better, so it is advisable to use both methods.

##  4 Experiments

In our experiments, we evaluate Bayesian Optimization using an encoder-based
warm-starts and compare it to the rank-based approach and lack of warm-start.
As mentioned in the introduction, this is one of the possible meta-tasks. We
focus on it due to its popularity in research [Bergstra et al., 2011,
Nicoletta Del Buono et al., 2020, Xiyuan Liu et al., 2020, Victoria and
Maragatham, 2021]. We use BO because, in the early stages of the process, it
has a problem with cold start. To mitigate that, one can provide a set of
starting points to speed up the algorithm’s convergence.

###  4.1 Data

We select three sets of datasets for the experiment purposes:  • [ (1)]  (1) A
classification tasks proposed in [Iwata and Kumagai, 2020] from OpenML
[Vanschoren et al., 2013] for training of the encoders.  (2) A subset of UCI
repository [Markelle Kelly, 2021] chosen by the authors of D2V [Jomaa et al.,
2021] for evaluation of warm-start on medium-sized data.  (3) A medical
dataset MIMIC-IV obtained from physionet.org for evaluation of warm-start
methods on small datasets. The latter set of datasets has 12 different tasks
which share all their features (both meanings and distributions).

Details. The set of datasets from OpenML are binary classification tasks
created from regression ones when a classification target is obtained as an
indicator of whether the regression target is above or below the average. The
datasets are preprocessed and standardized, so the distribution of ones and
zeros is equal. A complete list of these datasets is in Appendix A.

The next set of data is supplied from the GitHub repository of D2V. Some of
the datasets are multilabel classification tasks with a number of classes
varying from 2 to 10. However, we decide to cast classes so that there are
only 2 classes in every task. The reason for such modification is that the
tasks on which the encoders are trained are binary only. We also remove
datasets with more than 10000 observations or 20 features. Appendix B shows a
complete list of used datasets.

The last set of data is prepared on the metaMIMIC datasets [Woźnica et al.,
2023]. It contains fully preprocessed data from [Johnson et al., 2023], with
12 binary classification tasks, each describing the presence of a specific
disease. As tasks originally contain more than 100 features, selecting the 10
most important for each task provides the ability to measure the similarity of
tasks by simply calculating the fraction of shared most significant features.
During our study, we would like to consider several small datasets, as it is
the preferred data type for liltab. Due to that fact, we split each metaMIMIC
task into multiple smaller tasks. For every small task used in the warm-start
evaluation, we specify 4 observations as the training set and 29 observations
as the validation set. Both training and validation sets are sampled equally
with respect to the target value. With this methodology, we yielded over 5000
small tasks in total. MetaMIMIC has not yet been used in this domain,
especially as the source of small subsets of data. However, it is a promising
type of set of tasks for the dataset similarity as each task has some common
features with the others.

###  4.2 Training of encoders

We train both encoders on OpenML datasets. The details of the training process
are listed in Appendix C. Full reproduction of the training encoders, as well
as, the following experiments, can be done with code on GitHub:
https://github.com/azoz01/rethinking_encoder_warmstart.

###  4.3 Optimized Algorithms

During our research on the warm-start, we focus only on searching for the best
HP for the elastic net [Pedregosa et al., 2011] and XGBoost [Chen and
Guestrin, 2016] classifiers. We use these algorithms because the first is
simpler, and the other is more complex. Also, they are characterized by
different tunability [Probst et al., 2019]. Moreover, the elastic net is
proven to perform well on small and medium-sized datasets that we use in
experiments. On the other hand, XGBoost is one of the state-of-the-art models
on tabular data [Grinsztajn et al., 2022].

(a) The representations obtained with liltab.

(b) The representations obtained with D2V.

Figure 1: T-SNE visualization of three encoded datasets’ representations.
Datasets diabetes_diagnosed-10 and hypertensive_diagnosed-10 share 7 features
while hypotension_diagnosed-10 shares 3 features with the others. Each point
represents an encoded subset of data sampled from one of the datasets. This
figure presents that both models encode subsets of the same dataset close to
each other; thus, encoders fulfill Requirement 1.

###  4.4 Setup of experiments

For the experiment’s purpose, we generate a random set of the hyperparameters
of the elastic net and XGBoost, each consisting of 100 configurations. For
XGBoost, we use it only on UCI datasets, as metaMIMIC tasks would be too small
for this model. The definition of the ranges and distributions of the HP are
shown in Appendix D. Then, we train and evaluate the model with every
configuration of HP on each task. The evaluations on the HP configurations are
used in two ways: for encoders, hyperparameter configurations that were the
best on the $N$ closest in terms of Euclidean distance representations of
datasets are chosen as warm-start points. Warm-start points for the rank-based
method are chosen from the best on-average configurations from the net.

During the experiment, we compare methods of warm-start in terms of providing
the most gain in Bayesian Optimization iteration by iteration. The experiment
was performed on data from UCI and metaMIMIC. As a metric to measure gain, we
use Average Distance To Maximum (ADTM) [Yogatama and Mann, 2014] for scaled
ROC-AUC score. We use cross-validation with respect to the selection of
validation datasets. For metaMIMIC and UCI data, we use 4 folds (3 tasks in
validation fold) and 5 folds (20 tasks in validation fold), respectively. Each
optimization consists of 30 iterations with 10 warm-start iterations. The
experiments were run on a single CPU Intel Core i9-11950H with 2.6GHz clock
rate, 8 cores and 16 threads. We have 16 GB of RAM available.

(a) The representation obtained with liltab.

(b) The representation obtained with D2V.

Figure 2: T-SNE visualization of four encoded datasets’ representations.
Datasets are extracted from UCI and do not share any variables. Each point
represents an encoded subset of data sampled from one of the datasets. As in
the previous Figure, presented representations prove that models meet
Requirement 1.

##  5 Results

This section presents two different assessment methods of tabular data
representation. The first of them is the similarity of latent vectors
corresponding to similar datasets, i.e., meeting Requirement 1. The second one
is the metrics obtained from warm-start in HP optimization.

(a) Plots of ADTM of UCI data HPO with warm-start comparison. Lower score
means the performance of a method is closer to the best score obtained on a
specific dataset on average.

(b) Critical distance plot for logistic regression in 10th iteration. The
position on the scale denotes test statistic in Friedman test. Methods that
are connected with horizontal line are statistically indistinguishable.

(c) Critical distance plot for logistic regression in 30th iteration. The
position on the scale denotes test statistic in Friedman test. Methods that
are connected with horizontal line are statistically indistinguishable.

Figure 3: Results for UCI datasets, elastic net. Each fold in ADTM plots
present results for specific split of datasets to train and validation subset.
Critical distance plots show that the difference in performance between
methods is statistically irrelevant.

###  5.1 Encodings similarity of datasets

There is no explicit definition of dataset similarity and the best method to
create artificial datasets with fixed similarity is to sample a subset of
tasks (in terms of both observations as well as features and targets).
However, we take advantage of using the metaMIMIC dataset, as all 12 tasks in
this dataset share some of the features. Thus, the similarity of the two tasks
can be measured as the number of features they share. Here, we present results
that prove both encoder-based representations meet Requirement 1.

In Figure 1, we show the T-SNE representation of vectors transformed by liltab
and D2V. For this purpose, we choose three tasks: diabetes_diagnosis-10,
hypertensive_diagnosis-10 and hypotension_diagnosis-10. For both liltab and
D2V, there is a visible distinction between hypotension_diagnosis-10 and other
tasks. It is because diabetes_diagnosis-10 and hypertensive_diagnosed-10 share
7 out of 10 features, whereas hypotension_diagnosis-10 has 3 common features
with the remaining tasks.

In Figure 2, we show the T-SNE representation of UCI datasets. Both encoders
perform similarly, clustering 3 out of 4 datasets very well. In this case,
however, explaining why one dataset does not create a cluster is not as
straightforward as in metaMIMIC. Potentially, this dataset may share some
features with the others. However, such sharing is no longer related to the
semantics of the features but their distributions and relationships.

(a) Plots of ADTM of UCI data HPO with warm-start comparison. Lower score
means the performance of a method is closer to the best score obtained on a
specific dataset on average.

(b) Critical distance plot for XGBoost in 10th iteration. The position on the
scale denotes test statistic in Friedman test. Methods that are connected with
horizontal line are statistically indistinguishable.

(c) Critical distance plot for XGBoost in 30th iteration. The position on the
scale denotes test statistic in Friedman test. Methods that are connected with
horizontal line are statistically indistinguishable.

Figure 4: Results for UCI, XGBoost. Each fold in ADTM plots present results
for specific split of datasets to train and validation subset. Critical
distance plots show that the difference in performance between methods is
statistically irrelevant.

###  5.2 Transferability of HP

In this experiment, we show that Requirement 1 is not sufficient alone to
create representations that lead to the transferability of hyperparameters.
Meeting only this requirement might reflect some intrinsic complexity of
datasets, but this is not informative for the selection of a proper starting
point for HP optimization.

We present trajectories of ADTM for each fold with averages and confidence
intervals for each considered method. We also show the statistical
significance of our results aggregated across all folds on critical-distance
plots.

In Figure 3, we present the results for the elastic net on UCI datasets. Here,
right after the warm-start iteration, no method performs better than the one
without it. However, based on Figures 3(b) and 3(c), no method outperforms
others with statistical significance. Similar results are achieved in the case
of XGBoost on these datasets – no method occurs to be better than the others,
which can be observed in Figure 4.

In the case of the metaMIMIC datasets, there is a significant difference in
the results of the specific methods. Based on Figure 5(a), the rank-based
baseline is superior to the rest of the method. It is also shown in Figures
5(b) and 5(c): the rank method is significantly better than other methods.
Moreover, both encoders are indistinguishable in terms of performance on
$10$-th and $30$-th iteration.

(a) Plots of ADTM of metaMIMIC data HPO with warm-start comparison. Lower
score means the performance of a method is closer to the best score obtained
on a specific dataset on average.

(b) Critical distance plot for logistic regression in 10th iteration. The
position on the scale denotes test statistic in Friedman test. Methods that
are connected with horizontal line are statistically indistinguishable.

(c) Critical distance plot for logistic regression in 30th iteration. The
position on the scale denotes test statistic in Friedman test. Methods that
are connected with horizontal line are statistically indistinguishable.

Figure 5: Results for metaMIMC, elastic net. Each fold in ADTM plots present
results for specific split of datasets to train and validation subset. Here,
rank method results in the lowest average distance to maximum ROC-AUC score
for specific datasets. Critical distance plots show that the rank method
perform significantly better than other methods, while encoder-based methods
are indistinguishable in terms of performance.

##  6 Conclusion

In our study, we propose a new method to extract meaningful data
representation from the model proposed in [Iwata and Kumagai, 2020]. We
evaluate its performance and compare it to D2V, an established dataset encoder
[Zhu et al., 2023]. The experiments show that they have comparable results on
most warm-start tasks. In the visualizations, they present that obtained
dataset representations capture inter-datasets diversification and intra-
dataset similarity.

In the experiments, we compare the performance of both encoders to baselines.
The results show that the encoders perform comparable or worse than baselines
in the warm-start task, which can be seen on critical distance plots.

We use two significantly different collections of data for the experiment’s
purpose, challenging encoders to work on data that is not similar to the
training set of datasets. However, the data from metaMIMIC is out of
distribution for D2V, it has comparable results to liltab. Similarly, the UCI
datasets are adequate for D2V, but it has no better performance than liltab in
this case.

In our study, we show that using encoders has no significant gain over simpler
methods in the task of choosing warm-start points in BO. Moreover, our results
suggest that evaluated methods based on a meta-learning approach are not
necessarily better than simple random initialization.

The results suggest that meeting Requirement 1 does not imply transferability
HP between similar tasks. It may be caused by the fact that the encoders that
we use during our studies are optimized to create universal representations,
not connected to any specific meta-task. While this idea is intuitive,
creating encoders that impose a strict connection between their similarity and
HP transferability may be more advisable.

We also want to focus on the proposed rank-based method, which performs
comparably or better than the encoder. While simple, this method has some
advantages over more sophisticated methods, such as encoders. Firstly, in many
applications, some configurations are empirically proven to be the best in
most cases because of the low tunability of an algorithm [Probst et al.,
2019]. The rank-based method chooses the decent configurations for most tasks
that may outperform encoder-based methods when a new dataset differs from
already seen ones. Secondly, taking some good, yet independent of the current
tasks’ starting points may be considered an exploration over exploitation
solution. Methods that rely on the similarity of tasks, in cases when the
similarity does not imply HP transferability, choose bad starting points,
resulting in computing overhead in Bayesian Optimization.

##  7 Limitation and Broader Impact

Here, we want to emphasize some of the limitations of our work. Additionally,
we summarize the potential future directions in this field.

We only used two encoders because few can work on heterogeneous data.
Moreover, most of them have no public implementation. Considering the
benchmark provided in [Zhu et al., 2023], we choose the two best models that
work on varying feature and target spaces.

To perform experiments, we used three significantly different sets of
datasets. It is more than in most works in this field; however, it is
necessary to reproduce these experiments on a more extensive and versatile
collection of datasets to ensure that our results are not data-dependent.
Concerning reproducibility, one of the sets cannot be published (metaMIMIC).

We believe that creating new heuristics and encoders is the next direction in
this study. The results shown in this paper suggest that this field needs to
be explored more. In particular, there is a great need to create encoders or
heuristic-based methods that perform in all tasks better than random sampling
starting points.

Social impact statement

After careful reflection, we have determined that this work presents no
notable negative impacts to society or the environment.

## References

  * Iwata and Kumagai [2020] Tomoharu Iwata and Atsutoshi Kumagai.  Meta-learning from Tasks with Heterogeneous Attribute Spaces.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 6053–6063, 2020.
  * Vanschoren [2019] Joaquin Vanschoren.  Meta-learning.  In Frank Hutter, Lars Kotthoff, and Joaquin Vanschoren, editors, _Automated Machine Learning: Methods, Systems, Challenges_ , pages 35–61. Springer International Publishing, 2019.  doi: 10.1007/978-3-030-05318-5_2.
  * Jomaa et al. [2021] Hadi S Jomaa, Lars Schmidt-Thieme, and Josif Grabocka.  Dataset2Vec: Learning Dataset Meta-Features.  _Data Mining and Knowledge Discovery_ , 35:964–985, 2021.  doi: https://doi.org/10.1007/s10618-021-00737-9.
  * Davenport and Kalakota [2019] Thomas Davenport and Ravi Kalakota.  The potential for artificial intelligence in healthcare.  _Future Healthcare Journal_ , 6(2):94–98, 2019.  doi: 10.7861/futurehosp.6-2-94.
  * Alanazi [2022] Abdullah Alanazi.  Using machine learning for healthcare challenges and opportunities.  _Informatics in Medicine Unlocked_ , 30:100924, 2022.  doi: 10.1016/j.imu.2022.100924.
  * Rivolli et al. [2022] Adriano Rivolli, Luís P. F. Garcia, Carlos Soares, Joaquin Vanschoren, and André C. P. L. F. de Carvalho.  Meta-features for meta-learning.  _Knowledge-Based Systems_ , 240:108101, 2022.  doi: 10.1016/j.knosys.2021.108101.
  * Snoek et al. [2012] Jasper Snoek, Hugo Larochelle, and Ryan P Adams.  Practical Bayesian Optimization of Machine Learning Algorithms.  _Advances in Neural Information Processing Systems_ , 25, 2012.
  * Wistuba et al. [2015a] Martin Wistuba, Nicolas Schilling, and Lars Schmidt-Thieme.  Learning Hyperparameter Optimization Initializations.  In _2015 IEEE International Conference on Data Science and Advanced Analytics_. Institute of Electrical and Electronics Engineers Inc., 2015a.  doi: 10.1109/DSAA.2015.7344817.
  * Wistuba et al. [2015b] Martin Wistuba, Nicolas Schilling, and Lars Schmidt-Thieme.  Sequential Model-Free Hyperparameter Tuning.  In _2015 IEEE International Conference on Data Mining_ , pages 1033–1038, 2015b.  doi: 10.1109/ICDM.2015.20.
  * Wistuba et al. [2018] Martin Wistuba, Nicolas Schilling, and Lars Schmidt-Thieme.  Scalable Gaussian Process-Based Transfer Surrogates for Hyperparameter Optimization.  _Machine Learning_ , 107(1):43–78, 2018.  doi: 10.1007/s10994-017-5684-y.  Springer US.
  * Feurer et al. [2018] Matthias Feurer, Benjamin Letham, and Eytan Bakshy.  Scalable Meta-Learning for Bayesian Optimization.  2018\.  arXiv: 1802.02219.
  * Woźnica et al. [2023] Katarzyna Woźnica, Mateusz Grzyb, Zuzanna Trafas, and Przemys\law Biecek.  Consolidated learning: a domain-specific model-free optimization strategy with validation on metaMIMIC benchmarks.  _Machine Learning_ , pages 1–25, 2023.  doi: 10.1007/s10994-023-06359-0.  Publisher: Springer.
  * Vanschoren et al. [2013] Joaquin Vanschoren, Jan N. van Rijn, Bernd Bischl, and Luis Torgo.  OpenML: Networked Science in Machine Learning.  _SIGKDD Explorations_ , 15(2):49–60, 2013.  doi: 10.1145/2641190.2641198.
  * Markelle Kelly [2021] Kolby Nottingham Markelle Kelly, Rachel Longjohn.  UCI Machine Learning Repository, 2021.  URL https://archive.ics.uci.edu.
  * Song et al. [2015] Hyun Oh Song, Yu Xiang, Stefanie Jegelka, and Silvio Savarese.  Deep Metric Learning via Lifted Structured Feature Embedding.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 4004–4012, 2015.
  * Zaheer et al. [2017] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola.  Deep Sets.  In _Advances in Neural Information Processing Systems_ , volume 30, 2017.
  * Bergstra et al. [2011] James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl.  Algorithms for Hyper-Parameter Optimization.  In _Advances in Neural Information Processing Systems_ , volume 24, 2011.
  * Nicoletta Del Buono et al. [2020] Nicoletta Del Buono, Flavia Esposito, and Laura Selicato.  Methods for Hyperparameters Optimization in Learning Approaches: An Overview.  pages 100–112, 2020.  doi: 10.1007/978-3-030-64583-0_11.
  * Xiyuan Liu et al. [2020] Xiyuan Liu, XiYuan Liu, Jia Wu, Jia Wu, Senpeng Chen, and SenPeng Chen.  Efficient Hyperparameters optimization Through Model-based Reinforcement Learning and Meta-Learning.  _2020 IEEE 22nd International Conference on High Performance Computing and Communications; IEEE 18th International Conference on Smart City; IEEE 6th International Conference on Data Science and Systems (HPCC/SmartCity/DSS)_ , pages 1036–1041, 2020.  doi: 10.1109/hpcc-smartcity-dss50907.2020.00139.
  * Victoria and Maragatham [2021] A Helen Victoria and Ganesh Maragatham.  Automatic tuning of hyperparameters using Bayesian optimization.  _Evolving Systems_ , 12(1):217–223, 2021.  doi: 10.1007/s12530-020-09345-2.
  * Johnson et al. [2023] Alistair Johnson, Lucas Bulgarelli, Tom Pollard, Steven Horng, Leo Anthony Celi, and Roger Mark.  MIMIC-IV, 2023.
  * Pedregosa et al. [2011] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine Learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011.
  * Chen and Guestrin [2016] Tianqi Chen and Carlos Guestrin.  XGBoost: A Scalable Tree Boosting System.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 785–794, 2016.  doi: 10.1145/2939672.2939785.
  * Probst et al. [2019] Philipp Probst, Anne-Laure Boulesteix, and Bernd Bischl.  Tunability: Importance of Hyperparameters of Machine Learning Algorithms.  _The Journal of Machine Learning Research_ , 20(1):1934–1965, 2019.
  * Grinsztajn et al. [2022] Léo Grinsztajn, Edouard Oyallon, and Gaël Varoquaux.  Why do tree-based models still outperform deep learning on typical tabular data?  In _Advances in Neural Information Processing Systems_ , volume 35, pages 507–520, 2022.
  * Yogatama and Mann [2014] Dani Yogatama and Gideon Mann.  Efficient transfer learning method for automatic hyperparameter tuning.  In _Artificial intelligence and statistics_ , pages 1077–1085. PMLR, 2014.
  * Zhu et al. [2023] Max Zhu, Katarzyna Kobalczyk, Andrija Petrovic, Mladen Nikolic, Mihaela van der Schaar, Boris Delibasic, and Petro Lio.  Tabular few-shot generalization across heterogeneous feature spaces.  _arXiv preprint arXiv:2311.10051_ , 2023.

##  Appendix A Datasets from OpenML

Table 1: Datasets from OpenML used for encoder training. name | number of rows | number of columns
---|---|---
diggle_table_a1 | 48 | 5
fri_c4_100_25 | 100 | 26
chatfield_4 | 235 | 13
edm | 154 | 18
vineyard | 52 | 3
pollution | 60 | 16
fri_c3_100_10 | 100 | 11
visualizing_hamster | 73 | 6
pyrim | 74 | 28
visualizing_slope | 44 | 4
hip | 54 | 8
chscase_geyser1 | 222 | 3
rabe_148 | 66 | 6
bodyfat | 252 | 15
chscase_vine1 | 52 | 10
sleuth_ex1714 | 47 | 8
echocardiogram-uci | 132 | 8
hutsof99_child_witness | 42 | 16
fri_c3_250_10 | 250 | 11
chscase_funds | 185 | 2
qsartox | 16 | 24
qsabr2 | 15 | 10
sleep | 62 | 8
longley | 16 | 7
heart | 270 | 14
mu284 | 284 | 10
ICU | 200 | 20
qsbralks | 13 | 22
analcatdata_uktrainacc | 31 | 16
slump | 103 | 10
rabe_131 | 50 | 6
gascons | 27 | 5
rabe_265 | 51 | 7
rabe_166 | 40 | 2
baskball | 96 | 5
bolts | 40 | 7
chscase_demand | 27 | 11
qsfsr2 | 19 | 10
fri_c2_100_5 | 100 | 6
sleuth_ex1605 | 62 | 6
qsfsr1 | 20 | 10
visualizing_ethanol | 88 | 3
EgyptianSkulls | 150 | 5
transplant | 131 | 3
pwLinear | 200 | 11
diabetes_numeric | 43 | 3
autoPrice | 159 | 16
treepipit | 86 | 10
branin | 225 | 3
machine_cpu | 209 | 7
detroit | 13 | 14
fri_c2_250_25 | 250 | 26
qsbr_y2 | 25 | 10
humans_numeric | 75 | 15
visualizing_environmental | 111 | 4
USCrime | 47 | 14
rabe_176 | 70 | 4
fri_c0_250_5 | 250 | 6
rabe_266 | 120 | 3

##  Appendix B Datasets from UCI

Table 2: Datasets from UCI used for HPO warmstart evaluation. name | number of rows | number of columns
---|---|---
chess-krvkp | 3196 | 37
post-operative | 90 | 9
contrac | 1473 | 10
pittsburg-bridges-REL-L | 103 | 8
seeds | 210 | 8
ecoli | 336 | 8
acute-nephritis | 120 | 7
haberman-survival | 306 | 4
cylinder-bands | 512 | 36
pittsburg-bridges-TYPE | 105 | 8
ringnorm | 7400 | 21
tic-tac-toe | 958 | 10
led-display | 1000 | 8
hill-valley | 1212 | 101
waveform | 5000 | 22
credit-approval | 690 | 16
dermatology | 366 | 35
statlog-heart | 270 | 14
echocardiogram | 131 | 11
thyroid | 7200 | 22
planning | 182 | 13
spect | 265 | 23
musk-1 | 476 | 167
wall-following | 5456 | 25
ozone | 2536 | 73
hayes-roth | 160 | 4
breast-cancer-wisc-prog | 198 | 34
yeast | 1484 | 9
energy-y2 | 768 | 9
wine-quality-red | 1599 | 12
teaching | 151 | 6
musk-2 | 6598 | 167
zoo | 101 | 17
statlog-landsat | 6435 | 37
molec-biol-promoter | 106 | 58
heart-switzerland | 123 | 13
congressional-voting | 435 | 17
statlog-image | 2310 | 19
oocytes_trisopterus_states_5b | 912 | 33
oocytes_merluccius_nucleus_4d | 1022 | 42
wine-quality-white | 4898 | 12
mushroom | 8124 | 22
conn-bench-sonar-mines-rocks | 208 | 61
vertebral-column-3clases | 310 | 7
heart-hungarian | 294 | 13
low-res-spect | 531 | 101
semeion | 1593 | 257
iris | 150 | 5
ilpd-indian-liver | 583 | 10
optical | 5620 | 63
horse-colic | 368 | 26
waveform-noise | 5000 | 41
lenses | 24 | 5
glass | 214 | 10
mammographic | 961 | 6
twonorm | 7400 | 21
balance-scale | 625 | 5
abalone | 4177 | 9
pittsburg-bridges-T-OR-D | 102 | 8
bank | 4521 | 17
hepatitis | 155 | 20
breast-cancer-wisc | 699 | 10
ionosphere | 351 | 34
flags | 194 | 29
image-segmentation | 2310 | 19
breast-tissue | 106 | 10
monks-2 | 601 | 7
pittsburg-bridges-SPAN | 92 | 8
trains | 10 | 30
heart-cleveland | 303 | 14
spectf | 267 | 45
page-blocks | 5473 | 11
statlog-vehicle | 846 | 19
monks-3 | 554 | 7
blood | 748 | 5
oocytes_merluccius_states_2f | 1022 | 26
statlog-australian-credit | 690 | 15
energy-y1 | 768 | 9
heart-va | 200 | 13
steel-plates | 1941 | 28
breast-cancer | 286 | 10
lung-cancer | 32 | 57
wine | 178 | 14
spambase | 4601 | 58
oocytes_trisopterus_nucleus_2f | 912 | 26
pima | 768 | 9
annealing | 898 | 32
breast-cancer-wisc-diag | 569 | 31
lymphography | 148 | 19
car | 1728 | 7
cardiotocography-10clases | 2126 | 22
titanic | 2201 | 4
acute-inflammation | 120 | 7
pittsburg-bridges-MATERIAL | 106 | 8
monks-1 | 556 | 7
fertility | 100 | 10
balloons | 16 | 5
vertebral-column-2clases | 310 | 7
synthetic-control | 600 | 61
parkinsons | 195 | 23
statlog-german-credit | 1000 | 25
molec-biol-splice | 3190 | 61
cardiotocography-3clases | 2126 | 22

##  Appendix C Hyperparameters of encoders

Table 3: Hyperparameters of Liltab encoder Hyperparameter | Description | Value
---|---|---
num_epochs | Maximum number of epochs | 100000
learning_rate | Learning rate | 0.0001
weight_decay | Weight decay | 0
batch_size | Batch size | 37
early_stopping_epochs | Epochs without loss decrease before stop of training | 2500
hidden_representation_size | Dimensionality of output of intermediate networks | 32
n_hidden_layers | Number of hidden layers in intermediate networks | 3
hidden_size | Size of hidden layers in intermediate networks | 32
dropout_rate | Dropout rate in all intermediate networks | 0.1
Table 4: Hyperparameters of Dataset2Vec encoder Hyperparameter | Description | Value
---|---|---
gamma | Scaling factor in metric-based classification | 1
num_epochs | Maximum number of epochs | 100000
learning_rate | Learning rate | 0.001
weight_decay | Weight decay | 0.0001
batch_size | Batch size | 16
train_n_batches | Number of batches in one epoch | 100
early_stopping_epochs | Epochs without loss decrease before stop of training | 500
f_dense_hidden_size | Size of input/output layers in residual block in first layer | 16
f_res_hidden_size | Size of hidden layers in residual block in first layer | 16
f_res_n_hidden | Number of hidden layers in residual block in first layer | 3
f_dense_out_hidden_size | Output size of first layer | 16
f_block_repetitions | Number of residual blocks in first layer | 3
g_layers_sizes | Sizes of hidden layers in second layer | [32, 16, 8]
h_dense_hidden_size | Size of input/output layers in residual block in third layer | 32
h_res_hidden_size | Size of hidden layers in residual block in third layer | 32
h_res_n_hidden | Number of hidden layers in residual block in third layer | 3
h_dense_out_hidden_size | Output size of last layer | 16
h_block_repetitions | Number of residual blocks in third layer | 3

##  Appendix D Hyperparameters grids used for HPO

Table 5: Hyperparameters grid used for HPO of elastic net algorithm. Condition | Hyperparameter | Values’ range | Distribution
---|---|---|---
n/a | tol | [0.0001, 0.001] | loguniform
n/a | C | [0.0001, 10000] | loguniform
n/a | solver |  | [lbfgs, liblinear, newton-cg,
---
newton-cholesky, sag, saga]
categorical
solver = liblinear | intercept scaling | [0.001, 1] | uniform
penalty | [l1, l2] | categorical
| solver = liblinear
---
and penalty = l2
dual | [true, false] | categorical
solver = saga | penalty | [elasticnet, l1, l2, null] | categorical
l1 ratio | [0, 1] | uniform
|  solver $\neq$ saga
---
and solver $\neq$ liblinear
penalty | [l2, null] | categorical
Table 6: Hyperparameters grid used for HPO of XGBoost algorithm. Hyperparameter | Values’ range | Distribution
---|---|---
no of estimators | [1, 1000] | uniform (int)
learning rate | [0, 1] | uniform
booster | [gblinear, gbtree] | categorical
subsample | [0.5, 1] | uniform
max depth | [6, 15] | uniform (int)
min child weight | [2, 256] | uniform
colsample bytree | [0.2, 1] | uniform
colsample bylevel | [0.2, 1] | uniform


# Detecting dynamical changes in time series by using the Jensen Shannon
Divergence

D. M. Mateos1, L. Riveaud2, P. W. Lamberti2
1 Neuroscience and Mental Health Programme, Division of Neurology, Hospital
for Sick Children.
Institute of Medical Science and Department of Paediatrics, University of
Toronto, Toronto, Canada.
2 Facultad de Matemática, Astronomía, Física y Computación (FaMAF), CONICET,
and Universidad Nacional de Córdoba, Córdoba, Argentina
*<EMAIL_ADDRESS>
###### Abstract

Most of the time series in nature are a mixture of signals with deterministic
and random dynamics. Thus the distinction between these two characteristics
becomes important. Distinguishing between chaotic and aleatory signals is
difficult because they have a common wide-band power spectrum, a delta-like
autocorrelation function, and share other features as well. In general signals
are presented as continuous records and require to be discretized for being
analyzed. In this work we present different schemes for discretizing and for
detection of dynamical changes in time series. One of the main motivations is
to detect transition from chaotic regime to random regime. The tools used are
originated in Information Theory. The schemes proposed are applied to
simulated and real life signals, showing in all cases a high proficiency for
detecting changes in the dynamics of the associated time series.

##  1 Introduction

In many areas of science the dynamics of processes underlying the studied
phenomena can be represented by means of finite time series, that can be
interpreted as realizations of stochastic processes (SP) or as orbits of
deterministic processes. Although they are very different in nature, Wold
demonstrated that under certain restrictions there is a relationship between
them [1]. His famous theorem establishes that a stationary stochastic process
can be decomposed as a deterministic part (DP), which can be described
accurately by a linear combination of its own past, and a second part
represented as moving average component of finite order. A different situation
occurs when we have a non stationary time series. In these situations it is
not possible to separate the series in a deterministic and in a stochastic
part. Moreover, in cases in which the DP is chaotic, it is possible to find a
DP which produces a time series that could be very difficult to distinguish of
a SP [2, 3]. For these reasons, the issue of distinguishing between time
series produced by deterministic chaos from noise has led to the development
of different scheme of analysis. This problem has already been treated with
various techniques such as: complexity-entropy plane [4, 5, 6], Lyapunov
exponents [7, 2, 8] or applying neural networks [9], among others. Chaotic
signals always produce time series with strong physical structure, unlike
those originated in stochastic processes. These have a little structure
depending on their correlations factors [10]. The basic idea behind our
schemes is to use quantifiers of how close or far are two time series. The
quantifiers is the Jensen–Shannon Divergence (JSD), which was introduced as a
measure of similarity between two probability distributions function (PDF)
[11].
Quantifiers based on information theory, such as the JSD require to associate
to a given time series a probability distribution (PD). Thus the determination
of the most suitable PD is a fundamental problem, because of the PD and the
sample space are inextricably linked. Many methods have been proposed to deal
with this problem. Among others, we can mention binary symbolic dynamics [12],
Fourier analysis [13], frequency counting [14], wavelet transform [15] or
comparing the order of neighboring relative values [16]. In a similar
approach, we propose to explore a way of mapping continuous-state sequences
into discrete-state sequences via the so-called alphabetic mapping [17].
Here we present two methods for analyzing time series. The first one provides
a relative distance between two different time series. It consists in the
application of the JSD to a previously discretized signal, through alphabetic
mapping. The quantity so evaluated will be called the alphabetic
Jensen–Shannon divergence (aJSD). The second method detects the changes in the
dynamics of a time series, through an sliding window which run over the
previously mapped signal. This window is divided in two parts. To each part we
can assign the corresponding PD, which can be compared by using the JSD. In
both cases the change in the underlying dynamics of the time series can be
detected by analyzing the changes in the value of the JSD at the neighborhood
of the time when the change occurs. Both methods were tested by using well
known chaotic and random signals taken from the literature and biophysical and
mechanical real signals.
The structure of this work is the following: In section 2 we give a brief
introduction of the Jensen Shannon divergence, and we explained in detail the
alphabetic mapping. In section 3 we introduce the aJSD. In section 4 we give a
brief characterization of the chaotic maps and colored noises. We display the
analysis of the chaotic and noisy series in section 5. In section 6 we
illustrate how the alphabetic Jensen Shannon divergence can be used for real
data extracted from biological and mechanical systems. We finish the paper by
drawing up some conclusions.

##  2  Jensen-Shannon divergence for alphabetic mapping

###  2.1 Jensen-Shannon Divergence

The Jensen Shannon divergence (JSD) is a distance measure between probability
distributions introduced by C. Rao [18] and J. Lin [19] as a symmetric version
of the Kullback–Leibler divergence. I. Grosse and co-workers improved a
recursive segmentation algorithm for symbolic sequences based on the JSD [20].
Let $X$ be a discrete-states variable $x_{i}$, $i=1,2,\ldots,N$ and let
$P_{1}$ and $P_{2}$ be two probability distributions for $X$, which we denote
as $p_{i}^{(1)}=P_{1}(x_{i})$ and $p_{i}^{(2)}=P_{2}(x_{i})$, with $0\leq
p_{i}^{(k)}\leq 1$ and $\sum_{i=1}^{N}p_{i}^{(k)}=1$ for all $i=1,2,\ldots,n$
and $k=1,2$. If $\pi_{1}$ and $\pi_{2}$ denote the weights of $P_{1}$ and
$P_{2}$ respectively, with the restrictions $\pi_{1}+\pi_{2}=1$ and
$\pi_{1},\pi_{2}\geq 0$, the JSD is defined by

$D_{JS}(P_{1},P_{2})=H(\pi_{1}P_{1}+\pi_{2}P_{2})-(\pi_{1}H(P_{1})+\pi_{2}H(P_{2}))$
(1)

where

$H(P)=-C\sum_{i=1}^{N}p_{i}\;log(p_{i})$ (2)

is the Shannon entropy for the probability distribution $P$. We take
$C=1/log(2)$ so that the entropy (and also the JSD) is measured in bits.
It can be shown that the JSD is positive defined, symmetric and it is zero if
and only if $P_{1}=P_{2}$ [11]. Moreover JSD is always well defined, it is the
square of a metric [21] and it can be generalized to compare an arbitrary
number of distributions, in the following way: Let $P_{1}(x),...,P_{M}(x)$ be
a set of probability distributions and let $\pi_{1},...,\pi_{M}$ be a
collection of non negative numbers such that $\sum_{j=1}^{M}\pi_{j}=1$. Then
the JSD between the probability distributions $P_{j}(x)$ with $j=1,..,M$ is
defined by

$D_{JS}\left(P_{1},...,P_{M}\right)=H\left[\sum_{j=1}^{M}\pi_{j}P_{j}\right]-\sum_{j=1}^{M}\pi_{j}H\left[P_{j}\right].$
(3)

###  2.2 Alphabetic Mapping

In Nature it can be found discrete-state sequences, such as DNA sequences or
those generated by logical circuits, but most of the “ the real life” signals
are presented as continuous records. Thus to use the JSD first we have to
discrete the data 111 It is possible to calculate entropy for continuous
states signals but the estimation of a differential entropy from the data is
not an easy task [22]. . There are many ways to overcome this requirement.
Here we use the alphabetic mapping introduced by Yang [17].
For a given continuous series $\\{x_{t}\\}$, we can mapping this real-values
series in a binary series easily, depending the relative values between two
neighbor points $[x_{t};x_{t+1}]$ in the following way:

$s_{n}=\left\\{\begin{array}[]{l l}0,&\quad\text{if $x_{t}\leq x_{t+1}$}\\\
1,&\quad\text{if $x_{t}>x_{t+1}$}\end{array}\right.$

$S$ is a binary sequence $S=\\{s_{1},s_{2},...,s_{N-1}\\}$ where $s_{i}=0$ or
$1$. Consider two integers $d\geq 2$ and $\tau\geq 1$, let us then define a
trajectory in the $d$-dimensional space as

$\mathbf{Y}_{t}\rightarrow[s_{t-(d-1)\tau},...,s_{t}]\qquad t\geq(d-1)\tau$

The vector $\mathbf{Y}^{(d,\tau)}_{t}$ is called “alphabetic word´´, where $d$
is the “embedding dimension´´ (the number of bit taken to create the word) and
$\tau$ is called the ”delay”. Taken’s theorem gives conditions on $d$ and
$\tau$, such that ${\mathbf{Y}^{d,{\tau}}_{t}}$ preserves the dynamical
properties of the full dynamic system (e.g, reconstruction of strange
attractors) [23, 24].

By shifting one data point at time, the algorithm produces a collection of
bit-word $\\{\mathbf{Y}^{d,\tau}\\}$ over the whole series. Therefore, it is
plausible that the occurrence of these bit-word reflects the underlying
dynamics of the original time series. Different types of dynamics produce
different distributions on these $\mathbf{Y}^{d,\tau}$ series. We define
$w^{d,\tau}_{i}$, as the symbol corresponding to the word
$\mathbf{Y}^{d,\tau}_{i}$. From these we construct a new series
$W^{d,\tau}=\\{w^{d,\tau}\\}$ which quantify the original series. The number
of the different symbols (alphabet length) depends on the number of bit taken;
in this case is $2^{d}$.

To give an example of the mapping we can consider the series
$\\{3.5,4.8,3.6,2,1,4.1,3.7,8.5,10.4,8.9\\}$, which has a corresponding the
serie $S=\\{1,0,0,1,0,1,1,0\\}$. For the election of the parameter $d=4$ and
$\tau=1$, the first word to appear is $\mathbf{Y}^{4,1}_{1}=(1,0,0,1)$, the
second one is $\mathbf{Y}^{4,1}_{2}=(0,0,1,0)$, and the other three are
$\mathbf{Y}^{4,1}_{3}=(0,1,0,1)$, $\mathbf{Y}^{4,1}_{4}=(1,0,1,1)$ and
$\mathbf{Y}^{4,1}_{5}=(0,1,1,0)$.

Frequently it is necessary to process signals of two or more dimensions such
as bi-dimensional chaotic maps, polysomnography, EEG, etc. The components of
such signals are mostly coupled, given signal values depending not only on the
previous values but also on the values reached by the other signals. Therefore
by making a one-dimensional analysis we can lost some valuable information.
Taking the idea discussed above for one-dimensional signals (1D), we have
applied the same algorithm with a slight modification to analyze 2-dimensional
(2D) signals without loosing information.
For a given continuous 2D series $\mathbf{X}_{t}=(x(t),y(t))$, we assign a 1D
string, by the relative values between the two components vector at each time
$t$, $[x_{t};y_{t}]$ in the following way:

$S_{t}=\left\\{\begin{array}[]{lll}0&\mbox{si}&x_{t}\geq y_{t}\\\\[5.69054pt]
1&\mbox{si}&x_{t}<y_{t}\end{array}\right.$

##  3 JSD combined with the alphabetic mapping

Once discretized the signals we can approximate the PD by the frequency of
occurrence of the symbols $W^{d,\tau}$. From this PD, we can develop the
following analysis schemes:

  * •

Distance between two signals: Given two different sequences, for example, a
chaotic and a random sequences, we map each one by using the methodos
explained in section 2.2. This gives us two sets of symbols
$W^{(d,\tau)}=\\{w^{(d,\tau)}\\}$ and
$\widetilde{W}^{(d,\tau)}=\\{\widetilde{w}^{(d,\tau)}\\}$. We calculate the
frequency appearance of the symbols for both sequences $P^{W}=P(W^{d,\tau})$
and $P^{\widetilde{W}}=P(\widetilde{W}^{d,\tau})$. Finally, we calculate the
Jensen Shannon divergence between these two distributions
$D_{JS}^{W,\widetilde{W}}=D_{JS}(P^{W}\lvert P^{\widetilde{W}})$.

  * •

Sliding Window: We introduce an sliding window that moves over the symbolic
sequence corresponding to the original signal. The window has a width
$\Delta>0$ and the position $k$ (referring the position of the center of the
window over the sequence). For each position $k$, we can divide the window in
two sub windows, one to the left and the other to the right of the position
$k$. For both windows, we evaluate the frequency of occurrence of symbolic
patterns: to the right $P(W^{d,\tau}_{r})=P^{W}_{r}$ and to the left
$P(W^{d,\tau}_{l})=P^{w}_{l}$. Finally we evaluate the associate JSD,
$D_{JS}(k)=D_{JS}(P^{W}_{r}\lvert P^{W}_{l})$ as a function of the pointer
position $k$. The position where the maximum value of the JSD occurs,
$D_{JS_{max}}=max[D_{JS}(k)]$, it is interpreted as the place where a
significant change in the probability distribution patterns $W^{(d,\tau)}$,
has occurred. This change can be associated to a variation in the statistical
properties of the original signal. The only restriction in the election of
window’s width is that $\Delta$ must be greater than the number of possible
patterns generated by the alphabetic mapping $(\Delta>>2^{d})$ .

Two sequences with the same statistical properties, should lead to identical
probability distributions and therefore, the divergence between them should
take a very small value, close to zero but non zero. This fact is due to the
statistical fluctuations. The estimators for probability distributions
corresponding to sequences must be constructed and then the fluctuations of
this construction will yield to JSD values greater than zero. To address this
problem Grosse et al. [20] introduced a quantity called ”significance” which
allows to see if the values reached by the JSD are greater than the
statistical fluctuations. This amount depends on the sequence length and the
size of the alphabet of symbols used in the representation of the sequence. An
expression for the significance value has presented in the reference [20]. A
limitation of that expression is that it is valid only for an alphabet with no
more than five symbols. Therefore we must modify the criteria introduced by
Grosse et al, to identify values of the JSD that are genuinely above
statistical fluctuations. To do that, we proposed to calculate the aJSD on a
set of $Ns=10^{5}$ ensembles generated for each signal using the same
parameters but with different initial conditions. After quantifying the
signals we have $Ns$ differents sequences
$(W^{d,\tau}_{1},....,W^{d,\tau}_{Ns})$.

The first step is to calculate the aJSD between all sequences belonging to the
same group of segments, that we call “auto-aJSD´´,

$D_{JS}^{W}=D_{JS}(P^{W_{i}}~{}\lvert~{}P^{W_{j}}),~{}i,j=1,...,Ns,~{}for~{}i\neq
j$

Then we evaluate the average value $\mu^{W}=\langle D_{JS}^{W}\rangle$ over
all the sets of sequences with its respective standard deviation
($\sigma^{W}=\langle(D_{JS}(P^{W_{i}}\lvert
P^{W_{j}})-\mu^{W})^{2}\rangle^{1/2}$). The next step is the same, but using
the two group of signals to be compared. For all the aJSD values resulting
from all the signal we take the average $\mu^{W,\widetilde{W}}=\langle
D_{JS}(P^{W_{i}}\lvert P^{\widetilde{W_{i}})}\rangle$, with its respective
standard deviation ($\sigma^{W,\widetilde{W}}=\langle(D_{JS}(P^{W_{i}}\lvert
P^{\widetilde{W_{i}}})-\mu^{W,\widetilde{W}})^{2}\rangle^{1/2}$).

Finally two sequences are different (in the statistical sense) if the
inequality:

$\mu^{W,\widetilde{W}}-\sigma^{W,\widetilde{W}}\geq\text{max}[\mu^{W}+\sigma^{W},\mu^{\widetilde{W}}+\sigma^{\widetilde{W}}]$

is satisfied. If the values of the aJSD do not pass this criteria we say that
the two signals are not statistically distinguishable one from each other.

##  4 Characterization of chaotic maps and colored noises

To test our scheme of analysis we used sequences extracted from the
bibliography. We use 18 chaotic maps and 5 colored noises, that we describe
briefly.

###  4.1 Chaotic maps

We consider 18 chaotic maps which were taken from the reference [25]. They can
be grouped as follows:

  1. 1.

1D chaotic maps: also call non–inverted maps. They are dynamical sequences for
wich the image has more than one pre-image and in each interaction a loss of
information occurs, generating in this way a chaotic system [25].

     * •

the lineal congruential generator [26],

     * •

the gaussian map [27],

     * •

the logistic map [28],

     * •

the Pinchers map [29],

     * •

the Ricker’s population model [30],

     * •

the sine circle [31],

     * •

the sine map [32],

     * •

the Spencer map [33] and

     * •

the tent map [34].

The maps that have been used are :

     * •

the Hénon map [35],

     * •

the Lonzi map,

     * •

the delayed logistic map [36],

     * •

the tinkerbell map [37],

     * •

the dissipative standard map [38],

     * •

the Arnold’s cat map [39],

     * •

the chaotic web map [40],

     * •

the Chirikov standard map [41],

     * •

the gingerbreadman map [42] and

     * •

the Hénon area-preserving quadratic map [35].

###  4.2 Colored noises

Noises with a power spectrum that varies with frequency are called colored
noises. There are many types of noises, depending on the shape of the power
spectrum and the distribution of values. The noise power spectrum often varies
with frequency as $1/f^{\alpha}$ (some time called Hurst noise). White noise
correspond to $\alpha=0$. To generate all noises we used the algorithm
described in [4, 43]

##  5  Distinguishing between chaotic and random sequences

In this section we present the results obtained from evaluating the aJSD
between chaotic maps and colored noises. We apply the two methods described in
section 2. First we create an aJSD distance matrix between chaotic and noise
signals and between chaotic and chaotic maps. In the second place we merge two
different signals (i.e. one chaotic and one noisy) and through the aJSD
sliding window method it is possible to detect the signal changes from one
regime to another.

###  5.1  Distance matrix between sequences

For each type of process explained in the previous section we generate
$N_{s}=10^{6}$ time series of $L_{s}=10^{6}$ data points with identical
parameters [25] and a random initialation. We compute a distances matrix
between chaotic and colored noises, using the significance criterion for the
aJSD values as explained in section 3. For the discretization of the signals
we used the parameters $\tau=1$ and $8<d<12$. Fig. 1 displays the matrix for
the corresponding parameters $d=8$ and $\tau=1$. For different embedding
dimension, we obtained similar results.

In the case of the aJSD distances matrix corresponding to chaos-noise (Fig: 1)
we can observe that most of chaotic maps are distinguishable from the
different types of colorated noises. The number in the boxes represent the
values of the aJSD. Lower is this value, more similar are the sequences.

Only for the particular case of the lineal congruential generator map (LCG)
and white noise (WN), the aJSD value does not pass the significance criterion.
An interpretation of this is that the LCG map is an example of a random number
generator passing the Miller–Rabin test [44]. Therefore the distribution of
words $\\{W^{(8,1)}\\}$ corresponding to LCG and WN are similar.

Figure 1: aJSD distance matrix between 1D chaotic sequences and colored
signals. The values chosen of the parameters $d$ and tau are: $d=8$ and
$\tau=1$. It can be observed that there is a good discrimination between
chaotic sequences and noise. aJSD values vary according to the correlation
noise and the dynamics of chaotic series. The value 0.00 means that the aJSD
has not exceeded the significance criterion set out in section 3 so that the
two sequences are indistinguishable from each other. Similar results were
found for parameters $8<d<12$ and $\tau=1$, showing that the method is robust
respect the election parameter $d$.

The same evaluation have been done between chaotic maps. Figure 2 displays the
corresponding distance matrix. As it can be seen all the aJSD values are above
significance criterion, showing our method to be adequated to distinguish
between different types of chaos. A more detailed analysis shows a strong
relationship between the values of aJSD and the phase diagram of the chaotic
maps. For example, the logistic map and the tent map have similar phase
diagrams 222 All phase diagrams belong to chaotic time series used in this
work are shown in the Sprott book (Apendix [25]), and the value of the
distance between both sequences is very low. The same behaviour can be found
between the Pincher’s map and the Spencer map.

Figure 2: aJSD distance matrix between 1D chaotic sequences. The values chosen
for the parameters $d$ and $\tau$ are: $d=8$ and $\tau=1$. It can be observed
a good discrimination between the different chaotic sequences. The absolute
values of aJSD vary according to the similarity between the phase diagrams of
the chaotic sequences, giving information of the similarity about the dynamics
of such sequence. Similar dynamics give similar sequences. Similar results
were found for parameters $8<d<12$ and $\tau=1$

The same analysis was performed for 2D chaotic series. In this case we use the
two-dimensional assignment described in section 2.2 and where the parameters
$8<d<12$ and $\tau=1$ were chosen . Fig 3 shows the corresponding aJSD
distance matrix for the parameters $d=8$ and $\tau=1$. As observed for one-
dimensional chaotic maps, all distances over passed the significance criteria,
being all maps distinguishable one from each other. There exists an strong
correspondence between the values of the aJSD and the similarity (or
dissimilarity) of the topology of the phase diagrams of the maps. The aJSD
values decrease as the topology of the phase diagram tend to be similar, can
see it in the case of Hénon map and the Lonzi map. Conversely when two maps
are topographical different, such Hénon map preserver area and Chirikov map,
the aJSD increas. Different embedding dimensions changes the absolute value of
the aJSD but the relative value between the elements of the distance matrix
remain unaltered.

Let us recall that in our scheme the aJSD measure the distance between the PDF
associated with the set of words $\\{W^{(d,\tau)}\\}$ and
$\\{\widetilde{W}^{(d,\tau)}\\}$. As a consequence of the Taken’s theorem
[24], these sets of words are in correspondence with certain aspects of the
phase space of each original signal. Series which have a similar dynamics,
have similar phase space and PDF comparable giving low values of the aJSD.

Figure 3: aJSD distance matrix between 2D chaotic sequences. The values chosen
for the parameters $d$ and $\tau$ are: $d=8$ and $\tau=1$. The matrix shows
aJSD values is a good discrimination between the different chaotic sequences.
As happen in the 1D chaotic maps, when more similar are the phase space ,
lower is the value of aJSD. Similar results were found for parameters $8<d<12$
and $\tau=1$

###  5.2 Detection changes in time series

Here we use the proposed sliding window scheme for detecting changes in a
signal. For this purpose we used two different signals of equal length
$L_{x_{1}}=L_{x_{2}}=5.10^{4}$ symbols each one, which are merged in a single
sequence, where the signal $x_{1}$ is a chaotic one and $x_{2}$ is a random
one, or two different chaotic sequences. Examples of two combined normalized
sequences are plotted in Figures. 4 A and B. Figure 4 shows the results
obtained by applying the aJSD to four combinations of chaos-noise sequences.
The aJSD achieves its maximum value exactly at the merging point of the two
sequences, which is marked with a dotted vertical line. Similar results were
observed in the case of sequences generated by chaotic processes. In these
cases the aJSD value reaches several orders of magnitude higher than those
corresponding to a single stationary sequence.

Figure 4: (A) Example of a combined chaos-noise signal (Rikers maps - Blue
noise). The signals were merged at the meddle position (point=200). (B) idem A
but using a chaos-chaos combined signals (Sine Circle - Linear congruential
Map). (C) aJSD calculation using the running window method, applied to four
composed chaotic and random sequences. The values chosen for the parameters
$d$ , $\tau$ and $\Delta$ are: $d=8$, $\tau=1$ and $\Delta=40000$, . (D) Idem
C applied to four composed different chaotic signals. In both plots, the
maximum aJSD value is reached at the point where the two different sequences
were merged (dotted vertical line)

##  6 Application of the aJSD to real data

Now we test our proposed schemes for a group of real world signals. The first
one is a set of ECG signals on which we calculated the distance matrix of aJSD
between groups of patients with different cardiac pathologies. In the second
application we use our methods to detected the alignment of the axis of an
electric motor.

###  6.1 Distinction between groups of patients with heart diseases

In this study case the signals correspond to the time intervals between
heartbeats (BBI) on 15 patients. The patients were grouped into 3 set of 5
patients each one. The first group consists of healthy persons with normal
sinus rhythm (NSR); the second set contains patients suffering from congestive
heart failure (CHF), and the third one is composed of patients suffering
atrial fibrillation (AF). These set of data are freely available at:
www.physionet.org/challenge/chaos/.
Each series represents a record of 24 hours (approximately 100,000 intervals).
The analysed records do not have any previous filter. Figure 5 plots the aJSD
distance matrix among the three sets, for the parameters ( $d=8$ and
$\tau=1$). The aJSD can discriminate between the members of the control group
and the pathological groups, being greater the distance between the members of
the pathological groups and the members of the control group. The frequency
time between heartbeats for the patients belonging to the control group,
remains virtually unchanged, while for patients with pathologies the
interbeats time is altered. The difference between heartbeats is captured
through the alphabetical mapping, giving different patterns distributions of
$W^{d,\tau}$ for each group.

Figure 5:  aJSD distance matrix for three groups of 5 patients each one: the
normal sinus rhythm (NSR) group, the congestive heart failure (CHF) group and
the atrial fibrillation (AF) group. The aJSD can discriminate between this
tree groups.The parameter used in this matrix was $d=8$ and $\tau=1$. The same
results were obtained for higher embedding dimensions ($8\leq d\leq 12$) .

###  6.2 Misalignment detection of an electric motor

The second example corresponds to a record taken from the axis movement of an
electric motor. The first signal corresponds to the vibration measure using a
capacitive accelerometer, mounted along the $z$ axis on the motor bearing
(Figure 6A). The second one, is a signal taken from an optical incremental
rotary encoder drum, with about 145 pulses per revolution as shown in Figure (
6B). The third is a signal generated by a piezoelectric accelerometer mounted
in the same place that the capacitive signal, which is plotted in Figure (
6C). The last signal is the rate of the engine load and depicted in Figure
(6D) 333It is related to the frequency of rotation of the shaft by means of a
”Keyphaser” . The data were obtained with a sampling frequency of $25.10^{3}$
Hz without any preprocessing. For more technical details on the recording
setup see (reference [45])
Each signal has a length of $N=N_{a}+N_{m}=14.10^{4}$, in where the first
$N_{a}=7.10^{4}$ measurement values correspond to axis in the aligned
position and $N_{m}=7.10^{4}$ the in misaligned position. For all the signals
we used the following parameters: $d=8$, $\tau=1$ and $\Delta=4.10^{4}$. The
results are shown in Fig. 6. It should be noted that for all signals, the aJSD
maximum value is reached exactly when the shaft alignment state changes. That
point is identified by the dotted vertical line. In the case of the signals
from the vibration measurement by the piezoelectric accelerometers as
capacitive, the value of the maximum aJSD is smaller and more volatile than
the other two methods. This fact can be associated with the efficacy of the
measurement method. In all cases the different signals are clearly detectable
by the methods here proposed.

Figure 6: aJSD analysis using a sliding windows, corresponding to the four
different signals taken on the axis of an electric motor. (A) Signal of a
capacitive accelerometer mounted on the $z$ axis of the motor bearing. (B) An
optical signal incremental rotary encoder drum with 145 pulses per revolution.
(C) Signal generated by a piezoelectric accelerometer mounted on the same
place as A. (D) Signal engine load rate. All signals were analysed choosing
the parameters $d=8$, $\tau=1$, $\Delta=40000$.

##  7 Discussion

In paper we have introduced the notion of alphabetic Jensen Shannon Divergence
to measure a the distance between time series. The method can be used for one-
dimensional and two-dimensional signals. Using this measure, we have developed
two methods to distinguish between different types of signals. Particularly we
used them distinguish random from chaotic signals, and to compare different
types of chaotic signals between them. Through the aJSD distance matrix, we
were able to show that chaotic signals are clearly distinguishable from random
signals with diverse power spectrum. It is important to remark that chaotic
signals with different phase space structure are distinguishable from one
another, also. The values taken by the entries of the distance matrix, are in
correspond with the similarity (or dissimilarity of the corresponding phase
spaces). This fact was observed for one and two dimensional chaotic sequences.
Our methods were also applied to real world sequences. In particular some
corresponding to ECG records, from healthy and sick patients.
We also developed a procedure detect dynamical changes in a signal by using a
sliding window that moves on the resulting symbolic sequence after mapping the
original series by a alphabetical assignment, being the centre of the original
window moved a place per step. Then the aJSD is evaluated for the two
corresponding sub-sequences belonging to each half of the window. The maximum
value of the divergence corresponds to the point where occurs a change in the
probability distribution of the discretized signal. We tested the method using
two different signals merged. The aJSD allows to detect the point where the
two different signals were coupled in order to constitute signals with chaos-
noise and chaos-chaos parts. Finally we use this scheme to evaluate the
alignment condition of the axis of an electric motor.
Both methods have shown to be robust respect to the election of the parameters
in the probability distribution evaluation. Among the main advantage we can
highlight that both are easily implemented for the analysis of real world
signals.

## Acknowledgments

The authors thank Dr. Pierre Granjon from the GIPSA-lab, Grenoble, France, for
kindly providing them with the recording electric motor data and Consejo
Nacional de Investigaciones Científicas y Técnicas (CONICET) for financial
assistance. They also thank Dr. Jose Luis Perez–Velazquez for useful
discussions.

## References

  * [1] H. Wold.  A study in the analysis of stationary time series.  Almqvist & Wiksell, 1938.
  * [2] M. Cencini, M. Falcioni, E. Olbrich, H. Kantz, and A. Vulpiani.  Chaos or noise: Difficulties of a distinction.  Physical Review E, 62(1):427, 2000.
  * [3] H. Sakai and H. Tokumaru.  Autocorrelations of a certain chaos.  Acoustics, Speech and Signal Processing, IEEE Transactions on, 28(5):588–590, 1980.
  * [4] O. A. Rosso, H. A. Larrondo, M. T. Martin, A. Plastino, and M. A. Fuentes.  Distinguishing noise from chaos.  Physical review letters, 99(15):154102, 2007.
  * [5] L. Zunino, M. C. Soriano, and O. A Rosso.  Distinguishing chaotic and stochastic dynamics from time series by using a multiscale symbolic approach.  Physical Review E, 86(4):046210, 2012.
  * [6] F. Olivares, A. Plastino, and O. A. Rosso.  Contrasting chaos with noise via local versus global information quantifiers.  Physics Letters A, 376(19):1577–1583, 2012.
  * [7] W. A. Brock.  Distinguishing random and deterministic systems: Abridged version.  Journal of Economic Theory, 40(1):168–195, 1986.
  * [8] J. B. Gao, J. Hu, W. W. Tung, and Y. H. Cao.  Distinguishing chaos from noise by scale-dependent lyapunov exponent.  Physical Review E, 74(6):066204, 2006.
  * [9] J. B. Elsner.  Predicting time series using a neural network as a method of distinguishing chaos from noise.  Journal of Physics A: Mathematical and General, 25(4):843, 1992\.
  * [10] Z. Wu and N. E. Huang.  A study of the characteristics of white noise using the empirical mode decomposition method.  In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 460, pages 1597–1611. The Royal Society, 2004.
  * [11] T. M. Cover and J. A. Thomas.  Elements of information theory.  John Wiley & Sons, 2012.
  * [12] K. Mischaikow, M. Mrozek, J. Reiss, and A. Szymczak.  Construction of symbolic dynamics from experimental time series.  Physical Review Letters, 82(6):1144, 1999.
  * [13] G. E. Powell and I. C. Percival.  A spectral entropy method for distinguishing regular and irregular motion of hamiltonian systems.  Journal of Physics A, 12(11):2053, 1979.
  * [14] O. A. Rosso, H. Craig, and P. Moscato.  Shakespeare and other english renaissance authors as characterized by information theory complexity quantifiers.  Physica A: Statistical Mechanics and its Applications, 388(6):916–926, 2009.
  * [15] O. A. Rosso, S. Blanco, J. Yordanova, V. Kolev, A. Figliola, M. Schürmann, and E. Başar.  Wavelet entropy: a new tool for analysis of short duration brain electrical signals.  Journal of neuroscience methods, 105(1):65–75, 2001.
  * [16] C. Bandt and B. Pompe.  Permutation entropy: A natural complexity measure for time series.  Physical Review Letters, 88:174102, 2002.
  * [17] A. C. Yang, S. Hseu, H. Yien, A. L. Goldberger, and C. Peng.  Linguistic analysis of the human heartbeat using frequency and rank order statistics.  Physical review letters, 90(10):108103, 2003.
  * [18] C. R. Rao.  Differential metrics in probability spaces.  Differential geometry in statistical inference, 10:217–240, 1987\.
  * [19] J. Lin.  Divergence measures based on the shannon entropy.  Information Theory, IEEE Transactions on, 37(1):145–151, 1991.
  * [20] I. Grosse, P. Bernaola-Galván, P. Carpena, R. Román-Roldán, J. Oliver, and H. E. Stanley.  Analysis of symbolic sequences using the jensen-shannon divergence.  Physical Review E, 65(4):041905, 2002.
  * [21] D. M. Endres and J. E. Schindelin.  A new metric for probability distributions.  Information Theory, IEEE Transactions on, 49(7):1858–1860, 2003\.
  * [22] J. Beirlant, E. J Dudewicz, L. Györfi, and E. C. Van der Meulen.  Nonparametric entropy estimation: An overview.  International Journal of Mathematical and Statistical Sciences, 6(1):17–39, 1997.
  * [23] J. C. Robinson.  Dimensions, embeddings, and attractors, volume 186.  Cambridge University Press, 2011.
  * [24] F. Takens.  Detecting strange attractors in turbulence.  In Dynamical systems and turbulence, Warwick 1980, pages 366–381. Springer, 1981.
  * [25] J. C. Sprott.  Chaos and time-series analysis, volume 69.  Oxford University Press Oxford, 2003.
  * [26] D. E. Knuth.  Sorting and Searching, vol. 3.  1973\.
  * [27] W. H. Steeb and M. A. Van Wy.  Chaos and fractals: algorithms and computations.  Wissenschaftsverlag, 1992.
  * [28] R. M. May.  Simple mathematical models with very complicated dynamics.  Nature, 261(5560):459–467, 1976.
  * [29] A. Potapov and M. K. Ali.  Robust chaos in neural networks.  Physics Letters A, 277(6):310–322, 2000.
  * [30] W. E. Ricker.  Stock and recruitment.  Journal of the Fisheries Board of Canada, 11(5):559–623, 1954.
  * [31] V. I. Arnol’d.  Small denominators. i. mapping the circle onto itself.  Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, 25(1):21–86, 1961.
  * [32] S. H. Strogatz.  Nonlinear dynamics and chaos: with applications to physics, biology and chemistry.  Perseus publishing, 2001.
  * [33] Shaw R. and Z. Naturforsch.  Strange attractors, chaotic behavior, and information flow.  A, 36(80), 1981.
  * [34] M. W. Hirsch, S. Smale, and R. L. Devaney.  Differential equations, dynamical systems, and an introduction to chaos, volume 60.  Academic press, 2004.
  * [35] M. Hénon.  A two-dimensional mapping with a strange attractor.  Communications in Mathematical Physics, 50(1):69–77, 1976.
  * [36] D. G. Aronson, M. A. Chory, G. R. Hall, and R. P. McGehee.  Bifurcations from an invariant circle for two-parameter families of maps of the plane: a computer-assisted study.  Communications in Mathematical Physics, 83(3):303–354, 1982.
  * [37] H. E. Nusse.  Dynamics: numerical explorations, volume 101.  Springer, 1998.
  * [38] G. Schmidt and B. W. Wang.  Dissipative standard map.  Physical Review A, 32(5):2994, 1985.
  * [39] V. I. Arnolʹd and A. Avez.  Ergodic problems of classical mechanics, volume 9.  Benjamin, 1968.
  * [40] A. A. Chernikov, R. Z. Sagdeev, and G. M. Zaslavskii.  Chaos-how regular can it be?  Physics Today, 41:27–35, 1988.
  * [41] B. V. Chirikov.  A universal instability of many-dimensional oscillator systems.  Physics reports, 52(5):263–379, 1979.
  * [42] R. L. Devaney.  A piecewise linear model for the zones of instability of an area-preserving map.  Physica D: Nonlinear Phenomena, 10(3):387–393, 1984.
  * [43] Larrondo H.A.  program: noisef k.m, http://www.mathworks.com/matlabcentral/fileexchange/35381, 2012.
  * [44] M. O. Rabin.  Probabilistic algorithm for testing primality.  Journal of number theory, 12(1):128–138, 1980.
  * [45] M. Blödt, P. Granjon, B. Raison, and J. Regnier.  Mechanical fault detection in induction motor drives through stator current monitoring-theory and application examples.  Fault Detection, pages 451–488, 2010.