# A Statistically-Based Approach to Feedforward Neural Network Model Selection

Andrew McInerney***Department of Mathematics and Statistics, University of
Limerick<EMAIL_ADDRESS>Kevin Burke†††Department of Mathematics and
Statistics, University of Limerick<EMAIL_ADDRESS>
###### Abstract

Feedforward neural networks (FNNs) can be viewed as non-linear regression
models, where covariates enter the model through a combination of weighted
summations and non-linear functions. Although these models have some
similarities to the models typically used in statistical modelling, the
majority of neural network research has been conducted outside of the field of
statistics. This has resulted in a lack of statistically-based methodology,
and, in particular, there has been little emphasis on model parsimony.
Determining the input layer structure is analogous to variable selection,
while the structure for the hidden layer relates to model complexity. In
practice, neural network model selection is often carried out by comparing
models using out-of-sample performance. However, in contrast, the construction
of an associated likelihood function opens the door to information-criteria-
based variable and architecture selection. A novel model selection method,
which performs both input- and hidden-node selection, is proposed using the
Bayesian information criterion (BIC) for FNNs. The choice of BIC over out-of-
sample performance as the model selection objective function leads to an
increased probability of recovering the true model, while parsimoniously
achieving favourable out-of-sample performance. Simulation studies are used to
evaluate and justify the proposed method, and applications on real data are
investigated.

Keywords. Neural networks; Model selection; Variable selection; Information
criteria.

##  1 Introduction

Neural networks are a popular class of machine-learning models, which pervade
modern society through their use in many artificial-intelligence-based systems
(LeCun et al.,, 2015). Their success can be attributed to their predictive
performance in an array of complex problems (Abiodun et al.,, 2018). Recently,
neural networks have been used to perform tasks such as natural language
processing (Goldberg,, 2016), anomaly detection (Pang et al.,, 2021), and
image recognition (Voulodimos et al.,, 2018). Feedforward neural networks
(FNNs), which are a particular type of neural network, can be viewed as non-
linear regression models, and have some similarities to the models
traditionally used in statistical modelling (e.g., covariates enter the model
through a weighted summation, and the estimation of the weights for an FNN is
equivalent to the calculation of a vector-valued statistic) (Ripley,, 1994;
White,, 1989). There was a strong interest in neural networks within the
statistics community in the late 1980s and early 1990s as evidenced in the
comprehensive reviews by White, (1989), Ripley, (1993), and Cheng and
Titterington, (1994). Despite this early interest, however, these models have
still been historically less popular among statisticians, and the majority of
neural network research has been conducted outside of the field of statistics
(Breiman,, 2001; Hooker and Mentch,, 2021). Given this, there is a general
lack of statistically-based methods, such as model and variable selection,
which focus on developing parsimonious models.

Typically, the primary focus when implementing a neural network centres on
model predictivity (rather than parsimony); the models are viewed as ‘black-
boxes’ whose complexity is not of great concern (Efron,, 2020). It is perhaps
not surprising, therefore, that there is a tendency for neural networks to be
highly over-parameterised, miscalibrated, and unstable (Sun et al.,, 2022).
Nevertheless, FNNs can capture more complex covariate effects than is typical
within popular (linear/additive) statistical models. Consequently, there has
been renewed interest in merging statistical models and neural networks, for
example, in the context of flexible distributional regression (Rügamer et
al.,, 2020, 2021; Stasinopoulos et al.,, 2017) and mixed modelling (Tran et
al.,, 2020). However, statistically-based model selection procedures are
required to increase the utility of the FNN within the statistician’s toolbox.

Traditional statistical modelling is concerned with developing parsimonious
models, as it is crucial for the efficient estimation of covariate effects and
significance testing (Efron,, 2020). Indeed, model selection (which includes
variable selection) is one of the fundamental problems of statistical
modelling (Fisher,, 1922). It involves choosing the “best” model, from a range
of candidate models, by trading pure data fit against model complexity
(Anderson and Burnham,, 2004). As such, there has been a substantial amount of
research on model and variable selection (Miller,, 2002). As noted by Heinze
et al., (2018), typical approaches include significance testing combined with
forward selection or backward elimination (or a combination thereof);
information criteria such as AIC or BIC (Akaike,, 1973; Schwarz,, 1978;
Anderson and Burnham,, 2004); and penalised likelihood such as LASSO
(Tibshirani,, 1996; Fan and Lv,, 2010).

In machine learning, due to the focus on model predictivity, relatively less
emphasis is placed on finding a model that strikes a balance between
complexity and fit. Looking at FNNs in particular, the number of hidden nodes
is usually treated as a tunable hyperparameter (Bishop et al.,, 1995; Pontes
et al.,, 2016). Input-node selection is not as common, as the usual consensus
when fitting FNNs appears to be similar to the early opinion of Breiman,
(2001): “the more predictor variables, the more information”. However, there
are some approaches in this direction, and a survey of variable selection
techniques in machine learning can be found in Chandrashekar and Sahin,
(2014). In more complicated deep learning models, neural architecture searches
are sometimes used, which look to automate the architecture search procedure,
for example, through the use of evolutionary algorithms or reinforcement
learning (Elsken et al.,, 2019; Wistuba et al.,, 2019). However, when the
aforementioned approaches are implemented, the optimal model is usually
determined based on its predictive performance, such as out-of-sample mean
squared error, which can be calculated on a validation data set. Unlike an
information criteria, out-of-sample performance does not directly take account
of model complexity.

When framing an FNN statistically, there are several motivating reasons for a
model selection procedure that aims to obtain a parsimonious model. For
example, the estimation of parameters in a larger-than-required model results
in a loss in model efficiency, which, in turn, leads to less precise
estimates. Input-node selection, which is often ignored in the context of
neural networks, can provide the practitioner with insights on the importance
of covariates. Furthermore, eliminating irrelevant covariates can result in
cheaper models by reducing potential costs associated with data collection
(e.g. financial, time, energy). In this paper, we take a statistical view of
neural network selection by assuming an underlying (normal) error
distribution. Doing so enables us to construct a likelihood function, and,
hence, carry out information-criteria-based model selection, such as the BIC
(Schwarz,, 1978), naturally encapsulating the parsimony in the context of a
neural network. More specifically, we propose an algorithm that alternates
between selecting the hidden layer complexity and the inputs with the
objective of minimizing the BIC. We have found, in practice, that this leads
to more parsimonious neural network models than the more usual approach of
minimizing out-of-sample error, while also not compromising the out-of-sample
performance itself.

The remainder of this paper is structured as follows. In Section 2, we
introduce the FNN model while linking it to a normal likelihood function.
Section 3 motivates and details the proposed model selection procedure.
Simulation studies to investigate the performance of the proposed method, and
to compare it to other approaches, are given in Section 4. In Section 5, we
apply our method to real-data examples. Finally, we conclude in Section 6 with
a discussion.

##  2 Feedforward Neural Network

Let $y=(y_{1},y_{2},\dotsc,y_{n})\in\mathbb{R}^{n}$ be the response variable
of interest for a regression-based problem, where $n$ represents the number of
observations. For the $i$th observation, $i=1,\dotsc,n$, let
$x_{i}=(x_{1i},x_{2i},\dotsc,x_{pi})$ be a vector of $p$ covariates—the inputs
to the neural network model. We assume a model of the form

$y_{i}=g(x_{i})+\varepsilon_{i}$

where $\varepsilon_{i}$ is a random error that we assume has a
$N(0,\sigma^{2})$ distribution, and $g(\cdot)$ is a neural network,

$g(x_{i})=\gamma_{0}+\sum_{k=1}^{q}\gamma_{k}\phi\left(\sum_{j=0}^{p}\omega_{jk}x_{ji}\right).$
(2.1)

As we aim to frame FNNs as an alternative to other statistical non-linear
regression models (i.e., used on small-to-medium sized tabular data sets
relative to the much larger data sets seen more broadly in machine learning),
and due to the universal approximation theorem (Cybenko,, 1989; Hornik et
al.,, 1989), we are restricting our attention to FNNs with a single-hidden
layer. The parameters in Equation 2.1 are as follows: $\omega_{0k}$, the bias
term associated with the $k$th hidden node; $\omega_{jk}$, the weight that
connects the $j$th input node to the $k$th hidden node; $\gamma_{0}$, the bias
term associated with the output node; and $\gamma_{k}$, the weight that
connects the $k$th hidden node to the output node. The function $\phi(\cdot)$
is the activation function for the hidden layer, which is often a logistic
function. The number of parameters in the neural network is given by

$K=(p+2)q+1.$

A diagram of a neural network architecture with $p$ input nodes and $q$ hidden
nodes is shown in Figure 1. In the diagram, $x_{0}=1$, $h_{0}=1$, and

$h_{k}=\phi\left(\sum_{j=0}^{p}\omega_{jk}x_{ji}\right).$

Figure 1:  Neural network architecture with $p$ input nodes and $q$ hidden
nodes.

Given our assumption that $\varepsilon_{i}\sim N(0,\sigma^{2})$, we then make
use of the likelihood function

$\ell(\theta)=-\frac{n}{2}\log(2\pi\sigma^{2})-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(y_{i}-g(x_{i}))^{2},$

where
$\theta=(\omega_{01},\ldots,\omega_{p1},\ldots,\omega_{0q},\ldots,\omega_{pq},\gamma_{0},\ldots,\gamma_{q},\sigma^{2})^{T}$.
We maximise this likelihood to obtain $\hat{\theta}$ but note that the
estimates of the neural network parameters do not depend on the value of
$\sigma^{2}$, i.e., the residual sum of squares,
$\sum_{i=1}^{n}(y_{i}-g(x_{i}))^{2}$, can be estimated to obtain the neural
network parameters. This is useful since standard neural network software
(that minimises the residual sum of squares) such as nnet (Ripley and
Venables,, 2022) can be used used to optimise the neural network followed by
the estimation of $\sigma^{2}$ in a separate step.

The calculation of a likelihood function allows for the use of information
criteria when selecting a given model, and in particular, the Bayesian
information criterion (BIC) (Schwarz,, 1978),

$\text{BIC}=-2\ell(\hat{\theta})+\log(n)(K+1),$

where we have $K+1$ parameters, i.e., the $K$ neural network parameters plus
the variance parameter, $\sigma^{2}$. An attractive property of the BIC is
that it is “dimension-consistent”, i.e., the probability of selecting the
“true” model approaches one as sample size increases (Anderson and Burnham,,
2004)

##  3 Model Selection

To begin model selection, a set of candidate models must be considered. For
the input layer, we can have up to $p_{\text{max}}$ inputs, where
$p_{\text{max}}$ is the maximum number of covariates being considered, and
this is often the total number of covariates available in the data under
study. The input layer can contain any combination of these $p_{\text{max}}$
inputs. For the hidden layer, we must specify a $q_{\text{max}}$ value, which
is the maximum number of hidden nodes to be considered; this controls the
maximum level of complexity of the candidate models. We can then have between
one and $q_{\text{max}}$ nodes in the hidden layer. From a neural network
selection perspective, we aim to select a subset of $p\leq p_{\text{max}}$
covariates to enter the input layer and to build a hidden layer of $q\leq
q_{\text{max}}$ nodes to adapt to the required complexity. To carry out these
selections, we suggest a statistically-motivated procedure based on minimising
the BIC, since it directly penalises complexity and is known to be selection
consistent, i.e., BIC minimisation converges to the true model asymptotically.
In contrast, and more usually in machine learning applications, one could
consider predictive performance, for example, the out-of-sample mean squared
error. We will also consider this approach but find that it leads to
significantly more complex models than the use of BIC while only marginally
improving predictive performance.

###  3.1 Proposed Approach

We propose a stepwise procedure that starts with a hidden-node selection phase
followed by an input-node selection phase. (We find that this ordering leads
to improved model selection.) This is, in turn, followed by a fine-tuning
phase that alternates between the hidden and input layers for further
improvements. A schematic of this procedure is shown in Figure 2, and is
described at a high level in the following paragraphs.

Figure 2:  Model selection schematic. Nodes coloured grey are being considered
in current phase. Nodes coloured gold represent optimal nodes in that phase to
be brought forward to the next phase.

The procedure is initialised with the full model containing $p_{\text{max}}$
inputs and $q_{\text{max}}$ hidden nodes, and, as mentioned, starts with a
hidden-node selection phase. For each candidate model in this phase (i.e.,
models with $q\in\\{1,\ldots,q_{\text{max}}\\}$), the network optimiser is
supplied with $n_{\text{init}}$ random vectors of initial parameters, the log-
likelihood function is maximised at each of these vectors, and the overall
maximiser is found. The reason for supplying the neural network with different
vectors of initial parameters is due to the complex optimisation surface for
neural networks that may contain several local maxima. Thus, the use of a set
of initial vectors (rather than just one) aims to increase the chance of
finding the global maxima; of course, this cannot be guaranteed as is often
the case in more complex statistical models. Once all of the $q_{\text{max}}$
candidate models have been fitted, the hidden-node selection phase is
concluded by selecting the one whose hidden structure (i.e., number of nodes)
minimises the BIC.

Once the number of hidden nodes has been found, the hidden layer is fixed, and
model selection switches focus to the input layer; at this point, there are
$p_{\text{max}}$ inputs (i.e., covariates). For the input-node selection
phase, each input node is dropped in turn, with the aim of finding an input
whose removal yields a lower BIC; as with the previous phase, random sets of
initial parameters are used for each candidate model in the underlying
likelihood optimisation. If the removal of a given input node does yield a
lower BIC value, then that input node is dropped from the model (and if two or
more inputs result in a lower BIC, the one yielding the lower BIC is removed).
This is repeated until no covariate, when removed from the model, results in a
lower BIC.

Both the hidden layer and covariate selection phases are backward elimination
procedures. Rather than stopping the algorithm after these two phases, we find
it fruitful to search for an improved model in a neighbourhood of the current
“best” model by carrying out some further fine tuning. This is done by
considering the addition or removal of one hidden node, then the further
addition or removal of one input node, and these two steps are repeated
alternately until no further adjustment decreases the BIC. This fine-tuning
stage is analogous to stepwise model selection with backward and forward
steps. Note that one could apply this alternating stepwise procedure from the
offset, but we have found it to be significantly more computationally
efficient to focus first on the hidden and input layers (separately and in
that order) before moving to the stepwise phase.

The particular order of the model selection steps described above has been
chosen in order to have a higher probability in recovering the “true” model,
and to have a lower computational cost (see Section 4.1 for a detailed
simulation). Note that choosing the _set_ of input nodes requires a more
extensive search than choosing the _number_ of hidden nodes. There are more
candidate structures for the input layer as you can have any combination of
the nodes. Therefore, it is recommended to perform hidden-node selection
first, to eliminate any redundant hidden nodes and decrease the number of
parameters in the model before performing input-node selection.

##  4 Simulation Studies

In order to justify and evaluate the proposed model selection approach, a
number of simulation studies are used:

  * •

Simulation 1 (Section 4.1): In our first simulation study, we investigate the
effect of the ordering of the model selection steps to justify the procedure.
This includes the effect of performing input-node and hidden-node selection
phases first, the improvement of including a stepwise fine-tuning step, and
the performance of a procedure that only carries out iterative stepwise steps
(i.e., fine tuning from the offset).

  * •

Simulation 2 (Section 4.2): The second simulation study compares the
performance of using the BIC as the model selection objective function versus
using AIC or out-of-sample mean squared error (OOS).

  * •

Simulation 3 (Section 4.3): The third simulation study investigates the
performance of the proposed approach while varying the number of random
initialisations.

In each simulation study, the response is generated from an FNN with known
“true” architecture. The weights are generated so that there are three
important inputs, $x_{1},x_{2},x_{3}$, with non-zero weights, and ten
unimportant inputs, $x_{4},\ldots,x_{13}$, with zero weights. The “true”
hidden layer consists of $q=3$ hidden nodes, while we set our procedure to
consider a maximum of $q_{\text{max}}=10$ hidden nodes. In all simulation
studies, we vary the sample size $n\in\\{250,500,1000\\}$ and carry out 1,000
replicates. The metrics calculated to evaluate the performance of the model
selection approach are the average number of true zero weights _correctly_
dropped from the model (C), the probability of choosing the correct set of
inputs (PI), the probability of choosing the correct number of hidden nodes
(PH), and the probability of choosing the overall true model (PT). Our
proposed model selection approach is implemented in our publicly available R
package selectnn (McInerney and Burke,, 2022). The neural network function
used is nnet, which is available from the R package of the same name (Ripley
and Venables,, 2022).

###  4.1 Simulation 1: Model Selection Approach

This simulation study aims to justify the approach of the proposed model
selection procedure, i.e., a hidden-node phase, followed by an input-node
phase, followed by a fine-tuning phase; here, we label this approach as H-I-F.
Some other possibilities would be: to start with the input-node phase (I-H-F),
to stop the procedure without fine tuning (either H-I or I-H), or to only
carry out fine-tuning from the beginning (F). The objective function used for
model selection is BIC, and each approach has $n_{\text{init}}=5$ initial
vectors for the optimisation procedure. (The choice of objective function and
the effect of $n_{\text{init}}$ are investigated in Sections 4.2 and 4.3,
respectively.) The results of the simulation study are shown in Table 1.

Table 1:  Simulation 1: model selection metrics.

$n$ | Method | Time (s) | C(10) | PI | PH | PT
---|---|---|---|---|---|---
| H-I | 13.08 | 7.81 | 0.59 | 0.18 | 0.10
| I-H | 50.37 | 2.53 | 0.01 | 0.44 | 0.01
250 | H-I-F | 13.63 | 8.75 | 0.72 | 0.54 | 0.43
| I-H-F | 52.54 | 4.46 | 0.04 | 0.50 | 0.04
| F | 116.02 | 7.68 | 0.47 | 0.13 | 0.12
| H-I | 32.09 | 8.97 | 0.83 | 0.53 | 0.51
| I-H | 99.74 | 6.44 | 0.43 | 0.87 | 0.40
500 | H-I-F | 36.46 | 9.56 | 0.90 | 0.95 | 0.85
| I-H-F | 103.32 | 7.24 | 0.48 | 0.93 | 0.44
| F | 82.36 | 9.73 | 0.89 | 0.90 | 0.82
| H-I | 52.73 | 9.98 | 0.99 | 0.98 | 0.97
| I-H | 186.40 | 8.65 | 0.78 | 0.99 | 0.77
1000 | H-I-F | 52.66 | 9.99 | 0.99 | 1.00 | 0.99
| I-H-F | 188.81 | 8.80 | 0.78 | 1.00 | 0.77
| F | 168.70 | 9.89 | 0.97 | 0.99 | 0.96

  * Time (s), median time to completion in seconds (carried out on an Intel® Core™ i5-10210U Processor). Best values for a given sample size are highlighted in bold.

Comparing the methods without the fine-tuning stage, H-I has a greater
probability of recovering the true model versus I-H for all sample sizes.
However, when looking at layerwise selection, the probability of selecting the
correct structure is increased when that layer is selected in the second
phase, e.g., hidden-node selection is best when it comes second. This suggests
a relationship between the structure of the input and hidden layers (the
probability of correctly selecting the structure of one layer increases when
the other layer is more correctly specified). Therefore, H-I is likely better
than I-H due to input-node selection being a more difficult task than hidden-
node selection (determining the optimal _set_ of input nodes versus the
optimal _number_ of hidden nodes), and, hence, it is favourable to perform it
after hidden-node selection (given the number of hidden nodes is not
substantially larger than the number of input nodes).

The relationship between the structure of both layers can be handled by
incorporating a fine-tuning phase after both the H and I phases are completed.
Recall that the aim of fine tuning is to search for an improved solution in a
neighbourhood of the current solution, where both H and I steps are carried
out alternately (and include both backward and forward selections). Indeed, we
see that the addition of the fine-tuning phase improves on H-I in the smaller
sample sizes (in large part due to improved hidden-layer selection), but its
addition does not greatly improve on I-H. Moreover, the addition of fine
tuning only marginally adds to the computational expense. Overall, H-I-F is
significantly better than I-H-F both in terms of computational expense and
model selection.

One may also consider only carrying out fine-tuning steps from the offset,
which we denote by F. However, this does not perform as well as H-I-F at the
smallest sample size and is more computationally demanding.

From the above, the H-I-F approach is what we suggest as it leads to good
model selection performance while also being the most computationally
efficient approach.

###  4.2 Simulation 2: Model Selection Objective Function

This simulation study aims to determine the performance of using different
objective functions when carrying out model selection. Three objective
functions are investigated: BIC, AIC, and out-of-sample mean squared error
(OOS). The AIC approach is the same as the proposed approach in Section 3.1,
swapping BIC for AIC $=-2\ell(\hat{\theta})+2(K+1)$. The OOS approach follows
the same procedure, with the objective function replaced by out-of-sample mean
squared error, which is calculated on an additional validation data set that
is 20% the size of the training data set. As before, $n_{\text{init}}=5$
random initialisations are used. The results of the simulation study are shown
in Table 2.

Table 2:  Simulation 2: model selection metrics.

$n$ | Method | C(10) | PI | PH | PT
---|---|---|---|---|---
| AIC | 2.48 | 0.00 | 0.00 | 0.00
250 | BIC | 8.75 | 0.72 | 0.54 | 0.43
| OOS | 4.51 | 0.04 | 0.28 | 0.01
| AIC | 2.44 | 0.00 | 0.00 | 0.00
500 | BIC | 9.56 | 0.90 | 0.95 | 0.85
| OOS | 4.61 | 0.03 | 0.36 | 0.00
| AIC | 2.66 | 0.00 | 0.00 | 0.00
1000 | BIC | 9.99 | 0.99 | 1.00 | 0.99
| OOS | 5.29 | 0.02 | 0.46 | 0.00

  * Best values for a given sample size are highlighted in bold.

The results show that BIC far outperforms OOS and AIC in correctly identifying
the correct FNN architecture. Using OOS as the model selection objective
function never leads to correct neural network architecture being identified.
This is due to the inability of the OOS to correctly identify and remove the
unimportant covariates (C is always relatively low). Using AIC leads to even
worse performance, and this is likely due to the weaker penalty on model
complexity compared to BIC.

Given that different objective functions are being optimised, it is also of
interest to compare the approaches across a range of different criteria, e.g.,
compute the OOS and AIC values for models selected using BIC. Table 3 contains
the median AIC, BIC, and OOS values for the best model found by each approach
across all simulation replicates. In addition, the median out-of-sample mean
squared error evaluated on a test set is reported, which is computed on an
entirely new dataset (20% the size of the training set) that the OOS-
optimising procedure was not exposed to. Lastly, the median number of neural
network parameters, $K$, is also shown and note that the true value is $K=16$.
As expected, in general, minimising a given objective leads to the lowest
average value for that objective, e.g., a BIC-minimising procedure leads to
the lowest BIC value on average among the three procedures (AIC, BIC, OOS).
Interestingly however, the BIC procedure leads to very comparable OOS values
to the OOS procedure (and even beats the OOS procedure in the smallest samples
size); moreover, BIC-minimisation leads to the lowest OOS values on the test
data. This is particularly noteworthy since this is achieved using
approximately half as many parameters as the OOS-minimisation procedure. In
contrast, the models selected using AIC have worse out-of-sample performance
and significantly more parameters.

Table 3:  Simulation 2: model performance metrics.

$n$ | Method | AIC | BIC | OOS | OOS Test |  $K$ (16)
---|---|---|---|---|---|---
| AIC | 289.4 | 803.2 | 1.81 | 2.29 | 144
250 | BIC | 572.9 | 646.9 | 0.62 | 0.86 | 16
| OOS | 619.3 | 723.9 | 0.65 | 1.30 | 27
| AIC | 904.0 | 1507.7 | 0.98 | 1.03 | 144
500 | BIC | 1078.3 | 1152.9 | 0.51 | 0.53 | 16
| OOS | 1095.8 | 1267.2 | 0.50 | 0.57 | 37
| AIC | 2002.9 | 2695.5 | 0.70 | 0.76 | 141
1000 | BIC | 2142.3 | 2225.9 | 0.50 | 0.56 | 16
| OOS | 2157.6 | 2334.6 | 0.49 | 0.57 | 34

  * Lowest values for a given sample size are highlighted in bold.

###  4.3 Simulation 3: Number of Initialisations

Since FNNs have a complex optimisation surface (the log-likelihood function),
each model fit is supplied with $n_{\text{init}}$ random initial vectors with
the aim of avoiding local maxima. Of course, larger values of
$n_{\text{init}}$ improve the chances of finding the global maximum but
increase the computational expense. In the previous sections, we fixed
$n_{\text{init}}=5$, whereas, here, we vary it at
$n_{\text{init}}\in\\{1,5,10\\}$ using the proposed H-I-F BIC-minimisation
procedure. The results of the simulation study are shown in Table 4.

Table 4:  Simulation 3: model selection metrics.

$n$ | $n_{\text{init}}$ | Time (s) | C(10) | PI | PH | PT
---|---|---|---|---|---|---
| 1 | 1.90 | 8.02 | 0.52 | 0.17 | 0.06
250 | 5 | 13.63 | 8.75 | 0.72 | 0.54 | 0.43
| 10 | 32.30 | 8.87 | 0.72 | 0.89 | 0.68
| 1 | 9.37 | 7.35 | 0.44 | 0.28 | 0.12
500 | 5 | 36.46 | 9.56 | 0.90 | 0.95 | 0.85
| 10 | 70.69 | 9.78 | 0.95 | 0.99 | 0.93
| 1 | 12.98 | 8.93 | 0.62 | 0.64 | 0.38
1000 | 5 | 52.66 | 9.99 | 0.99 | 1.00 | 0.99
| 10 | 117.32 | 10.00 | 1.00 | 1.00 | 1.00

  * Time (s), median time to completion in seconds (carried out on an Intel® Core™ i5-10210U Processor). Best values for a given sample size are highlighted in bold.

We can clearly see the trade-off between better model selection and worse
computational efficiency as $n_{\text{init}}$ increases. We would certainly
recommend $n_{\text{init}}>1$ initial vectors since the results are poor for
$n_{\text{init}}=1$. Beyond this, the choice might be based on the
computational constraints in a given practical setting, but we note, in
particular, that larger values of $n_{\text{init}}$ are more important in
smaller sample sizes.

##  5 Application to Data

We consider two real-data applications for our proposed model selection
approach. For comparison purposes, the model found by our proposed model
selection procedure is compared to fitting an FNN with all inputs and the
maximum number of hidden nodes considered. Each data set has been randomly
split into a training set and test set with a 90%–10% split, respectively, and
all variables have been rescaled (based on the training data) to be in the
range $[0,1]$. The model selection procedure was implemented with
$n_{\text{init}}=10$ and $q_{\text{max}}=10$. For both models (selected and
full), we report the number of input nodes ($p$), the number of hidden nodes
($q$), the total number of parameters ($K$), the BIC, and the out-of-sample
mean squared error (OOS). The results are summarised in Table 5 and discussed
in more detail for each dataset in the Section 5.1 and 5.2.

###  5.1 Boston Housing Data

The Boston housing data set is available in James et al., (2022). It is a
modified version of a data set that originally came from a study by Harrison
Jr and Rubinfeld, (1978). It contains information relating to housing in 506
communities in the Boston area in 1970. The aim of the study was to examine
the relationship between twelve explanatory variables and the median house
price for each community. The explanatory variables for each community are the
per capita crime rate (crim), the proportion of residential land zoned for
lots over 25,000 sq.ft. (zn), the proportion of non-retail business acres per
town (indus), the average number of rooms per dwelling (rm), the weighted mean
of distances to five Boston employment centres (dis), the index of
accessibility to radial highways (rad), the pupil-to-teacher ratio (ptratio),
the proportion of owner-occupied units built prior to 1940 (age), the nitrogen
oxide concentration (parts per 10 million) (nox), the full-value property-tax
rate per $10,000 (tax), an indicator if the community bounds the Charles River
(chas), and the proportion of the population that are ‘lower status’ (lstat).
The response variable is the median value of owner-occupied homes ($1,000s).

The model selection procedure was implemented with $n_{\text{init}}=10$ and
$q_{\text{max}}=10$. The model found by the proposed approach has dropped
three covariates (chas, indus and rad), and found the optimal number of hidden
nodes to be five. It has a much lower BIC value than the full model
($\Delta\text{BIC}=187.8$), and an equivalent out-of-sample performance, while
having 85 fewer parameters.

###  5.2 Wine Quality Data Set

The wine quality data set is available on the UC Irvine Machine Learning
Repository, and was originally used by Cortez et al., (2009). The goal of the
data is to examine the relationship between wine quality and some
physicochemical tests. Due to the taste of red and white wine being quite
different, the data is split into two separate data sets: one for white wine
samples ($n$ = 4,898) and another for red wine samples ($n$ = 1,599). There
are eleven predictor variables in total: fixed acidity (g(tartaric
acid)/$\text{dm}^{3}$) (fixed), volatile acidity (g(acetic
acid)/$\text{dm}^{3}$) (vol), citric acid (g/$\text{dm}^{3}$) (cit), residual
sugar (g/$\text{dm}^{3}$) (sugar), chlorides (g(sodium
chloride)/$\text{dm}^{3}$) (chl), free sulfur dioxide (mg/$\text{dm}^{3}$)
(fsulfur), total sulfur dioxide (mg/$\text{dm}^{3}$) (tsulfur), density
(g/$\text{cm}^{3}$) (den), pH (pH), sulphates (g(potassium
sulphate)/$\text{dm}^{3}$) (sulphate), and alcohol (vol.%) (alco). The
response variable is the quality of the wine, which was graded on a scale from
zero (very bad) to ten (excellent) by sensory assessors.

For the red wine data, the model found by the proposed approach dropped four
covariates: fixed, sugar, fsulfur, and den. The optimal number hidden nodes is
two. Again, the model selection procedure found a model with a much lower BIC
than the full model ($\Delta\text{BIC}=298.5$). It has an equivalent out-of-
sample performance to the full model with 112 fewer parameters.

For the white wine data, the model selected dropped two covariates, chl and
alco, and selected $q=4$. Again, the model selected has lower BIC than the
full model ($\Delta\text{BIC}=234.1$)., and has 86 fewer parameters with an
equivalent mean squared error.

Table 5:  Application Results

| Boston | Red Wine | White Wine
---|---|---|---
| Selected | Full | Selected | Full | Selected | Full
$p$ | 9 | 12 | 7 | 11 | 9 | 11
$q$ | 5 | 10 | 2 | 10 | 4 | 10
$K$ | 56 | 141 | 19 | 131 | 45 | 131
BIC | -1145.8 | -958 | -1953.5 | -1655.0 | -6132.1 | -5898.0
OOS | 0.01 | 0.01 | 0.02 | 0.02 | 0.01 | 0.01

  * $p$, number of input nodes; $q$, number of hidden nodes; $K$, number of parameters

##  6 Discussion

FNNs have become very popular in recent years and have the potential to
capture more complex covariate effects than traditional statistical models.
However, model selection procedures are of the utmost importance in the
context of FNNs since their flexibility may increase the chance of over-
fitting; indeed, the principle of parsimony is very common throughout
statistical modelling more generally. Therefore, we have proposed a
statistically-motivated neural network selection procedure by assuming an
underlying (normal) error distribution that then permits BIC minimisation.
More specifically, our procedure involves a hidden-node selection phase,
followed by an input-node (covariate) selection phase, followed by a final
fine-tuning phase. We have made this procedure available in our selectnn
package in R (McInerney and Burke,, 2022).

Through extensive simulation studies, we have found that that (i) the order of
selection (input versus hidden layer) is important, with respect to the
probability of recovering the true model and the computational efficiency,
(ii) the addition of a fine-tuning stage provides a non-negligible improvement
while not significantly increasing the computational burden, (iii) using the
BIC is necessary to asymptotically converge to the true model, and (iv)
although the models selected using BIC have fewer parameters than those
selected using out-of-sample performance, they have comparable, and sometimes
improved, predictivity. We suggest that statistically-orientated model
selection approaches are necessary in the application of neural networks —
just as they are in the application of more traditional statistical models —
and we have demonstrated the favourable performance of our proposal.

We expect that randomisation and/or divide-and-conquer throughout the
selection phases would be required in more complex problems involving many
covariates and/or hidden layers. Nevertheless, neural networks are still
valuable in more traditional (smaller) statistical problems for which
procedures such as ours will lead to more insightful outputs. Furthermore, the
implementation of statistical approaches more broadly (such as uncertainty
quantification and hypothesis testing) in neural network modelling will be
crucial for the enhancement of these insights. This will be the direction of
our future work.

## Acknowledgement

This work has emanated from research conducted with the financial support of
Science Foundation Ireland under Grant Number 18/CRT/6049.

## References

  * Abiodun et al., (2018) Abiodun, O. I., Jantan, A., Omolara, A. E., Dada, K. V., Mohamed, N. A., and Arshad, H. (2018).  State-of-the-art in artificial neural network applications: A survey.  Heliyon, 4(11):e00938.
  * Akaike, (1973) Akaike, H. (1973).  Information theory and an extension of the maximum likelihood principle.  In Selected Papers of Hirotugu Akaike, pages 199–213. Springer.
  * Anderson and Burnham, (2004) Anderson, D. and Burnham, K. (2004).  Model selection and multi-model inference.  Second. NY: Springer-Verlag.
  * Bishop et al., (1995) Bishop, C. M. et al. (1995).  Neural networks for pattern recognition, chapter 9, pages 353–354.  Oxford university press.
  * Breiman, (2001) Breiman, L. (2001).  Statistical modeling: The two cultures (with comments and a rejoinder by the author).  Statistical science, 16(3):199–231.
  * Chandrashekar and Sahin, (2014) Chandrashekar, G. and Sahin, F. (2014).  A survey on feature selection methods.  Computers & Electrical Engineering, 40(1):16–28.
  * Cheng and Titterington, (1994) Cheng, B. and Titterington, D. M. (1994).  Neural networks: A review from a statistical perspective.  Statistical Science, pages 2–30.
  * Cortez et al., (2009) Cortez, P., Cerdeira, A., Almeida, F., Matos, T., and Reis, J. (2009).  Modeling wine preferences by data mining from physicochemical properties.  Decision Support Systems, 47(4):547–553.
  * Cybenko, (1989) Cybenko, G. (1989).  Approximation by superpositions of a sigmoidal function.  Mathematics of Control, Signals and Systems, 2(4):303–314.
  * Efron, (2020) Efron, B. (2020).  Prediction, estimation, and attribution.  International Statistical Review, 88:S28–S59.
  * Elsken et al., (2019) Elsken, T., Metzen, J. H., and Hutter, F. (2019).  Neural architecture search: A survey.  The Journal of Machine Learning Research, 20(1):1997–2017.
  * Fan and Lv, (2010) Fan, J. and Lv, J. (2010).  A selective overview of variable selection in high dimensional feature space.  Statistica Sinica, 20(1):101.
  * Fisher, (1922) Fisher, R. A. (1922).  On the mathematical foundations of theoretical statistics.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 222(594-604):309–368.
  * Goldberg, (2016) Goldberg, Y. (2016).  A primer on neural network models for natural language processing.  Journal of Artificial Intelligence Research, 57:345–420.
  * Harrison Jr and Rubinfeld, (1978) Harrison Jr, D. and Rubinfeld, D. L. (1978).  Hedonic housing prices and the demand for clean air.  Journal of Environmental Economics and Management, 5(1):81–102.
  * Heinze et al., (2018) Heinze, G., Wallisch, C., and Dunkler, D. (2018).  Variable selection–a review and recommendations for the practicing statistician.  Biometrical Journal, 60(3):431–449.
  * Hooker and Mentch, (2021) Hooker, G. and Mentch, L. (2021).  Bridging breiman’s brook: From algorithmic modeling to statistical learning.  Observational Studies, 7(1):107–125.
  * Hornik et al., (1989) Hornik, K., Stinchcombe, M., and White, H. (1989).  Multilayer feedforward networks are universal approximators.  Neural Networks, 2(5):359–366.
  * James et al., (2022) James, G., Witten, D., Hastie, T., and Tibshirani, R. (2022).  ISLR2: Introduction to Statistical Learning, Second Edition.  R package version 1.3-1.
  * LeCun et al., (2015) LeCun, Y., Bengio, Y., and Hinton, G. (2015).  Deep learning.  Nature, 521(7553):436–444.
  * McInerney and Burke, (2022) McInerney, A. and Burke, K. (2022).  selectnn: A Statistically-Based Approach to Neural Network Model Selection.  R package version 0.0.0.9000 (https://github.com/andrew-mcinerney/selectnn).
  * Miller, (2002) Miller, A. (2002).  Subset selection in regression.  chapman and hall/CRC.
  * Pang et al., (2021) Pang, G., Shen, C., Cao, L., and Hengel, A. V. D. (2021).  Deep learning for anomaly detection: A review.  ACM Computing Surveys (CSUR), 54(2):1–38.
  * Pontes et al., (2016) Pontes, F. J., Amorim, G., Balestrassi, P. P., Paiva, A., and Ferreira, J. R. (2016).  Design of experiments and focused grid search for neural network parameter optimization.  Neurocomputing, 186:22–34.
  * Ripley and Venables, (2022) Ripley, B. and Venables, W. (2022).  nnet: Feed-forward neural networks and multinomial log-linear models.  R package version, 7.3-17.
  * Ripley, (1993) Ripley, B. D. (1993).  Statistical aspects of neural networks.  In Nielsen, B. O. E., Jensen, J. L., and Kendall, W. S., editors, Networks and Chaos: Statistical and Probabilistic Aspects, pages 40–123. Chapman & Hall.
  * Ripley, (1994) Ripley, B. D. (1994).  Neural networks and related methods for classification.  Journal of the Royal Statistical Society: Series B (Methodological), 56(3):409–437.
  * Rügamer et al., (2021) Rügamer, D., Kolb, C., Fritz, C., Pfisterer, F., Kopper, P., Bischl, B., Shen, R., Bukas, C., Sousa, L. B. d. A. e., Thalmeier, D., Baumann, P., Kook, L., Klein, N., and Müller, C. L. (2021).  deepregression: a flexible neural network framework for semi-structured deep distributional regression.  arXiv preprint arXiv:2104.02705.
  * Rügamer et al., (2020) Rügamer, D., Kolb, C., and Klein, N. (2020).  Semi-structured deep distributional regression: Combining structured additive models and deep learning.  arXiv preprint arXiv:2002.05777.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  Annals of statistics, 6(2):461–464.
  * Stasinopoulos et al., (2017) Stasinopoulos, M. D., Rigby, R. A., Heller, G. Z., Voudouris, V., and De Bastiani, F. (2017).  Flexible regression and smoothing: using GAMLSS in R.  CRC Press.
  * Sun et al., (2022) Sun, Y., Song, Q., and Liang, F. (2022).  Learning sparse deep neural networks with a spike-and-slab prior.  Statistics & Probability Letters, 180:109246.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288.
  * Tran et al., (2020) Tran, M.-N., Nguyen, N., Nott, D., and Kohn, R. (2020).  Bayesian deep net glm and glmm.  Journal of Computational and Graphical Statistics, 29(1):97–113.
  * Voulodimos et al., (2018) Voulodimos, A., Doulamis, N., Doulamis, A., and Protopapadakis, E. (2018).  Deep learning for computer vision: A brief review.  Computational Intelligence and Neuroscience, 2018.
  * White, (1989) White, H. (1989).  Learning in artificial neural networks: A statistical perspective.  Neural Computation, 1(4):425–464.
  * Wistuba et al., (2019) Wistuba, M., Rawat, A., and Pedapati, T. (2019).  A survey on neural architecture search.  arXiv preprint arXiv:1905.01392.


# Revisiting the Toda-Brumer-Duff criterion for order-chaos transition
in dynamical systems

F. Sattin<EMAIL_ADDRESS>Consorzio RFX (CNR, ENEA, INFN, Università
di Padova, Acciaierie Venete SpA), Corso Stati Uniti 4, 35127 Padova, Italy
L. Salasnich<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia
“Galileo Galilei” and QTech, Università di Padova, via Marzolo 8, 35131
Padova, Italy  Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via
Marzolo 8, 35131 Padova, Italy  Istituto Nazionale di Ottica del Consiglio
Nazionale delle Ricerche, Unità di Sesto Fiorentino, via Carrara 1, 50019
Sesto Fiorentino (Firenze), Italy

###### Abstract

The Toda-Brumer-Duff (TBD) is an analytical criterion for estimating the local
exponential rate of divergence between nearby trajectories in dynamical
systems, and it is employed as a test for assessing the existence of chaos
therein. It is fairly simple, intuitive, and works well in several situations,
hence gained quite a wide popularity, yet it is known to be not rigorous since
predicts “false positives”, i.e., flags as chaotic systems that are instead
regular. We revisit here the TBD criterion in order to understand the causes
of its failures, and pinpoint that the problem with it is due to two reasons:
(a) the TBD criterion does not constrain the trajectories to lie on the same
energy hypersurface; (b) it does not distinguish between the divergence of
trajectories along or perpendicularly to the direction of the flow, the former
being irrelevant for assessing the presence of chaos. We show how both points
can be incorporated within the TBD framework, yielding an amended criterion
which, when applied to reference cases, interprets correctly the kind of
dynamics observed. However, since must rely upon the knowledge of actual
trajectories, it cannot be employed as a shortcut with respect to the
numerical integration of the equations of motion.

Keywords: Toda-Brumer-Duff criterion; Hamiltonian systems; chaos;
integrability; tangent dynamics

##  I Introduction

Dynamical systems may feature both regular or irregular chaotic motion
depending on control parameters. Quite often one such parameter is the energy
$E$: at low $E$ the system is usually regular, and a smooth or sharp
transition to irregular motion occurs by increasing $E$. Identifying the kind
of behavior of a dynamical system is obviously extremely important; it can be
done directly through the estimate of the Lyapunov exponents or the Poincaré
plots, which are the two most known and employed numerical techniques. In the
course of time, several other algorithms have been developed, including the
Smaller Alignment Index (SALI) [1], the analysis of the power spectrum of
geodesic divergence between nearby trajectories [2], the Fast Lyapunov
Indicator method [3], as well as several others. All these techniques require
the numerical integration of the equations of motion for several trajectories.
In the past, when computer resources were more scarce, there was a premium for
overriding this computationally expensive step and assessing the
regular/irregular behavior through the analytical inspection of the structure
of the equations of motion. Even nowadays, analytical approaches may be
valuable since provide insight into the mechanisms that generate chaos and the
paths that drive the order/chaos transition. One such algorithm, valid for
two-dimensional systems not too far from integrable ones, is provided by
Codaccioni, Doveil and Escande [4].

The first analytical approach was developed in 1970 by Okubo [5], within the
framework of fluid mechanics, but it is nowadays more widely known from the
works of Toda [6], Brumer and Duff [7], applied to generic Hamiltonian
dynamical systems. Within this work, we will be restricted to Hamiltonians of
the most common kind, with separated kinetic and potential terms (and with
unitary mass for simplicity)

$H={|{\bf p}|^{2}\over 2}+V({\bf q})\;,$ (1)

where the potential $V({\bf q})$ is a nonlinear function of the $N$
coordinates ${\bf q}=(q_{1},q_{2},...,q_{N})$ and ${\bf p}=(p_{1},...,p_{N})$
is the vector of linear momenta. Hamilton equations write

$\displaystyle{\dot{\bf q}}(t)$ $\displaystyle=$ $\displaystyle{\bf p}(t)\;,$
(2) $\displaystyle{\dot{\bf p}}(t)$ $\displaystyle=$
$\displaystyle-\partial_{{\bf q}}V({\bf q}(t))\;,$ (3)

where the dot means the derivative with respect to time $t$ and $\partial_{\bf
q}$ means the gradient with respect to ${\bf q}$.
The stability of the motion is investigated by looking at two nearby
trajectories $({\bf q}_{1},{\bf p}_{1})$ and $({\bf q}_{2},{\bf p}_{2})$ and
rewriting the Hamilton equations in terms of the displacements $\delta{\bf
q}={\bf q}_{1}-{\bf q}_{2}$, $\delta{\bf p}={\bf p}_{1}-{\bf p}_{2}$ as
follows

$\displaystyle\delta{\dot{\bf q}}(t)$ $\displaystyle=$
$\displaystyle\delta{\bf p}(t)$ (4) $\displaystyle\delta{\dot{\bf p}}(t)$
$\displaystyle=$ $\displaystyle-\left[\partial_{\bf q}V(\delta{\bf q}(t)+{\bf
q}_{2}(t))-\partial_{\bf q}V({\bf q}_{2}(t))\right]\;.$ (5)

Defining the $2N$-dimensional displacement vector ${\bf D}=(\delta{\bf
q},\delta{\bf p})$, the right-hand sides of the previous equations are
linearized around ${\bf D}=0$:

${\dot{\bf D}}(t)={\hat{\bf M}(t)}\ {\bf D}(t)\;,$ (6)

where

$\displaystyle\hat{\bf M}(t)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&\hat{\bf I}_{N}\\\ -\hat{\bf H}&0\end{pmatrix}$
(7) $\displaystyle(\hat{\bf I}_{N})_{ij}$ $\displaystyle=$
$\displaystyle\delta_{ij}\;,\quad(\hat{\bf H})_{ij}=\partial^{2}_{ij}V({\bf
q}(t))$

Notice that $\hat{\bf M}$ depends from time only implicitly, through the
coordinates ${\bf q}$’s which evolve by virtue of Eqns. (2, 3). The solution
of the coupled system (2,3,6) over long times yields the largest Lyapunov
exponent of the system, a positive value flagging the existence of chaos. As
it is well known, the estimate of Lyapunov exponents can only be done
numerically. Toda, Brumer and Duff (TBD) postulate that a scale separation
holds: the coordinates of the common motion $({\bf q},{\bf p})$ are considered
as frozen while evolving ${\bf D}(t)$. That is, the entries of the $2N\times
2N$ matrix $\hat{\bf M}$ are considered as constant while solving Eq. (6),
which decouples from Eqns (2,3) and has exponential solutions

${\bf D}(t)=\sum_{i=1}^{2N}\chi_{i}{\bf d}_{i}\exp(\mu_{i}t)$ (8)

where $\mu_{i}$ are the eigenvalues of the matrix $\hat{\bf M}$, ${\bf d}_{i}$
the corresponding eigenvectors and $\chi_{i}$ numerical coefficients depending
upon initial conditions. Now the coordinates ${\bf q}$’s appearing in the
entries of the matrix $\hat{\bf M}$ (and therefore also in $\mu_{i}$) play the
role of arbitrary parameters, have no necessary reference to the particle’s
trajectory. Positive real-valued eigenvalues $\mu_{i}$ imply local exponential
divergence between neighboring trajectories at the position $\bf q$. TBD argue
that the regions where exponential divergence exists are the seeds for the
global chaotic behavior of the system. The rationale is that orbits visiting
repeatedly regions where at least one $\mu_{i}>0$ become unstable and
ultimately chaotic.

The TBD criterion is extremely intuitive and simple to implement; furthermore,
it proved effective in predicting order-chaos transitions for several systems,
thus it became quite popular. Besides generic Hamiltonian systems [8, 9, 10],
it was applied to problems from field theory [11, 12, 13], atomic physics
[14], nuclear physics [15, 16], plasma physics [17], neutral fluid and MHD
dynamics [18, 19, 20]. It was soon realized, though, that it is not a rigorous
criterion. Benettin et al. [21] showed that it does not exist a clear
correspondence between the loci of instability predicted by the TBD criterion
and the actual exponential departure of numerically computed trajectories:
chaos is a global property of a system, not a local attribute. Furthermore,
the TBD criterion predicts “false positives”, i.e., regions of hyperbolic
instability in systems that can be rigorously proved to be regular, including
one-dimensional systems like the pendulum, or two-dimensional systems like the
anti-Hènon-Heiles system [21, 22, 23]. These difficulties are way more
remarkable given the streamlined derivation of the TBD criterion. It is
obvious that the problem arises when decoupling Eq. (6) from Eqns. (2,3), but
it is not easy to identify where it lies. The purpose of the present paper is
to provide an answer to this question, with the hope of producing an amended
criterion free from these difficulties, leading to improved predictions.

##  II Scrutiny of the TBD criterion

The realization that the TBD criterion was not exact led immediately to
attempts to improve it. The first proposal was put forth by Cerjan and
Reinhardt [24]. They argued that the linearization step leading to Eq. (6)
actually entails investigating the stability properties of the fixed point
$\dot{\bf D}=0$ at ${\bf D}=0$. If, besides this one, further stable fixed
points at finite values $\bf D$ do exist, any initially unstable perturbation
will be ultimately dragged into the attraction basin of one of these points as
it grows. The Cerjan-Reinhardt suggestion is clever but, basically, is a model
for the nonlinear saturation of an instability. It does not explain why the
initial stages of the instability do not produce observable effects, in the
first place.

In 1981, Kosloff and Rice [25] realized that the displacement vector ${\bf D}$
within the TBD criterion is totally arbitrary, whereas two conditions are to
be fulfilled: (i) Both trajectories $({\bf q}_{1}(t),{\bf p}_{1}(t))$ and
$({\bf q}_{2}(t),{\bf p}_{2}(t))$ lie on the constant energy hypersurface
$H=E$, thus ${\bf D}(t)$ must be tangent to the same hypersurface. (ii) The
vector ${\bf D}(t)$ may be decomposed as

${\bf D}(t)={\bf D}_{\parallel}(t)+{\bf D}_{\perp}(t)\;,$ (9)

where ${\bf D}_{\parallel}(t)$ is parallel to the flow: ${\bf
D}_{\parallel}(t)\parallel(\dot{\bf q},\dot{\bf p})$, and ${\bf D}_{\perp}(t)$
is perpendicular to it (still lying on the surface $H=E$). Only the
exponential growth of ${\bf D}_{\perp}(t)$ is relevant in establishing the
chaotic properties of trajectories: choosing two initial points along the flow
is equivalent to considering one and the same trajectory, just starting at
different times. Within the TBD criterion, there is no way for discriminating
whether exponential growth takes place along ${\bf D}_{\parallel}(t)$ or ${\bf
D}_{\perp}(t)$.

It is straightforward to realize that accounting for the two conditions (i)
and (ii) rule out all one-dimensional systems ($N=1$) from the consideration
of the TBD criterion, since it is not possible to write couples of isoenergy
trajectories having at the same time ${\bf D}_{\perp}(t)\neq{\bf 0}$.

In the following, we will consider the specific but important case of two-
dimensional systems ($N=2$). We will show that in that case, too, a correct
consideration of the two conditions (i, ii) leads to a modified TBD criterion.

###  II.1 Two-dimensional systems: $N=2$

From now on we will use the symbols $(x,y)$ rather than $q_{1},q_{2}$, in
adherence with customary usage and in order not to burden notations.

Let us define the following local orthonormal basis system in the phase-space:

$\displaystyle{\bf e}_{0}$ $\displaystyle=$
$\displaystyle\left(\partial_{x}V,\partial_{y}V,p_{x},p_{y}\right)Q^{-1}$ (10)
$\displaystyle{\bf e}_{1}$ $\displaystyle=$
$\displaystyle\left(p_{x},p_{y},-\partial_{x}V,-\partial_{y}V\right)Q^{-1}$
(11) $\displaystyle{\bf e}_{2}$ $\displaystyle=$
$\displaystyle\left(-p_{y},p_{x},-\partial_{y}V,\partial_{x}V\right)Q^{-1}$
(12) $\displaystyle{\bf e}_{3}$ $\displaystyle=$
$\displaystyle\left(\partial_{y}V,-\partial_{x}V,-p_{y},p_{x}\right)Q^{-1}$
(13)

where $Q=\sqrt{(\partial_{x}V)^{2}+(\partial_{y}V)^{2}+{\cal{K}}^{2}}$, and
${\cal{K}}=\sqrt{p_{x}^{2}+p_{y}^{2}}=\sqrt{2(E-V)}$. Energy surfaces are
defined by the condition

$E=H={p_{x}^{2}\over 2}+{p_{y}^{2}\over 2}+V(x,y)$ (14)

The direction perpendicular to the energy surface is defined by the phase-
space gradient of $H:(\partial_{q}H,\partial_{p}H)$; it is straightforward to
check that this direction is spanned by ${\bf e}_{0}$. The second vector is
aligned along the flow: ${\bf
e}_{1}\propto\left({\dot{x}},{\dot{y}},{\dot{p}}_{x},{\dot{p}}_{y}\right)=\left(p_{x},p_{y},-\partial_{x}V,-\partial_{y}V\right)$.
Finally, the two vectors ${\bf e}_{2},{\bf e}_{3}$ span the remaining two
directions perpendicular to the flow.

We write the displacement vector as a linear combination over the basis (13),
restricting to the hyperplane tangent to the energy surface: that is, no
component along ${\bf e}_{0}$ is allowed. Thus, we write ${\bf D}$ in Eq. (6)
to be of the form

${\bf D}(t)=\sum_{i=1}^{3}\alpha_{i}(t){\bf e}_{i}.$ (15)

Consistently with the logic of the TBD approach, we keep the basis ${\bf
e}_{i}$ as fixed and let the coefficients $\alpha_{i}$ alone vary in time. We
replace this expression into Eq. (6) which project onto the basis vectors
${\bf e}_{1},{\bf e}_{2},{\bf e}_{3}$. We introduce the array
$\bm{\alpha}=(\alpha_{1},\alpha_{2},\alpha_{3})$ and the $3\times 3$ matrix
$\hat{\bf W}$: ${W}_{ij}={\bf e}_{i}\cdot\hat{\bf M}\cdot{\bf e}_{j}$. Thus,
we started from the $4\times 4$ array equation (6) and arrived to the $3\times
3$ one

${\dot{\bm{\alpha}}}(t)=\hat{\bf W}\ {\bm{\alpha}}(t)$ (16)

We reiterate that the entries of the matrix $\hat{\bf W}$ depend on the phase-
space coordinates $(x,y,p_{x},p_{y})$ which must here be considered as
parameters. The solution of (16) with given initial condition
$\bm{\alpha}(t=0)=\bm{\alpha}_{0}$ is

$\bm{\alpha}(t)=\sum_{i=1}^{3}(\bm{\alpha}_{0}\cdot{\bf
w}^{(i)})\,\exp(\lambda^{(i)}\,t)\,{\bf w}^{(i)}$ (17)

The $\lambda^{(i)}$’s are the eigenvalues of $\hat{\bf W}$, roots of the
third-order polynomial

$det(\hat{\bf W}-\lambda\hat{\bf
I})=-\lambda^{3}+c_{2}\lambda^{2}+c_{1}\lambda+c_{0}=0$ (18)

with

$\displaystyle c_{2}$ $\displaystyle=$
$\displaystyle{p_{x}\partial_{x}V(\partial_{xx}V-1)+p_{y}\partial_{y}V(\partial_{yy}V-1)+\partial_{xy}V(p_{y}\partial_{x}V+p_{x}\partial_{y}V)\over
Q^{2}},$ $\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle-{\partial_{xx}V(p_{y}^{2}+\partial_{y}V^{2})+\partial_{yy}V(p_{x}^{2}+\partial_{x}V^{2})-2\partial_{xy}V(p_{x}p_{y}+\partial_{x}V\partial_{y}V)\over
Q^{3}},$ $\displaystyle c_{0}$ $\displaystyle=$
$\displaystyle-{p_{x}\partial_{x}V(\partial_{xy}V^{2}+\partial_{yy}V-\partial_{xx}V\partial_{yy}V)+p_{y}\partial_{y}V(\partial_{xy}V^{2}+\partial_{xx}V-\partial_{xx}V\partial_{yy}V)-\partial_{xy}V(p_{y}\partial_{x}V+p_{x}\partial_{y}V)\over
Q^{4}}$

where $Q$ was defined in correspondence of Eqns. (13). Finally, the ${\bf
w}^{(i)}$’s are the eigenvectors associated to the eigenvalues
$\lambda^{(i)}$’s. The eigenvalues are the roots of the third-order polynomial
with real-valued coefficients (18), thus we know that one of the $\lambda$’s
must be purely real (but it may be either negative or positive). We will call
it $\Lambda_{R}$. The other two eigenvalues are complex-valued, conjugated to
each other:
$\Lambda_{I}=\Lambda_{I,r}+i\Lambda_{I,i},\Lambda_{I}^{*}=\Lambda_{I,r}-i\Lambda_{I,i}$.
The three corresponding eigenvectors will be labeled as ${\bf w}_{R},{\bf
w}_{I},{\bf w}_{I}^{*}$. In conclusion, Eq. (17) becomes

$\displaystyle\bm{\alpha}(t)$ $\displaystyle=$
$\displaystyle(\bm{\alpha}_{0}\cdot{\bf w}_{R})\exp(\Lambda_{R}\,t){\bf
w}_{R}+2Re\left[(\bm{\alpha}_{0}\cdot{\bf w}_{I})\exp(\Lambda_{I}\,t){\bf
w}_{I}\right]$ $\displaystyle=$ $\displaystyle(\bm{\alpha}_{0}\cdot{\bf
w}_{R})\exp(\Lambda_{R}t){\bf
w}_{R}+2\exp(\Lambda_{I,r}\,t)\,\left\\{\cos(\Lambda_{I,i}\,t)Re\left[(\bm{\alpha}_{0}\cdot{\bf
w}_{I}){\bf
w}_{I}\right]-\sin(\Lambda_{I,i}\,t)Im\left[(\bm{\alpha}_{0}\cdot{\bf
w}_{I}){\bf w}_{I}\right]\right\\}$

In order to make further progress, we recall the second point raised by
Kosloff and Rice: The existence of real-valued positive eigenvalues is not a
sufficient condition to infer instability. We can envisage two potential
instabilities corresponding to either $\Lambda_{R}>0$ or $Re(\Lambda_{I})>0$.
The second case is not meaningful, though, since the presence of oscillating
terms in $(\Lambda_{I,i}\,t)$ to a large extent averages towards zero the
corresponding dilatation. The only meaningful case to consider, therefore, is
that when, for some triplet $(x,y,\phi)$, $\Lambda_{R}>0$. The term
proportional to $\exp(\Lambda_{R}\,t)$, therefore, dominates over the others
in Eq. (17) (but for particular choices of the initial conditions), thus the
whole sum collapses onto just this single term

$\bm{\alpha}(t)\simeq(\bm{\alpha}_{0}\cdot{\bf
w}_{R})\exp(\Lambda_{R}\,t)\,{\bf w}_{R}$ (20)

We insert this expression in Eq. (15):

${\bf D}(t)=\sum_{i=1}^{3}\alpha_{i}(t){\bf
e}_{i}\simeq\sum_{i=1}^{3}(\bm{\alpha}_{0}\cdot{\bf
w}_{R})\exp(\Lambda_{R}\,t)\,({\bf w}_{R})_{i}\,{\bf e}_{i}$ (21)

Since it may generate confusion, it is worthwhile recalling that $\bf D$,
${\bf e}_{i}$ are vectors in the physical phase-space, they have units.
Conversely, $\bm{\alpha}$, ${\bf w}_{R,I}$ are vectors in an abstract space:
they are pure numbers. In (21) the expression $\sum_{i=1}^{3}({\bf
w}_{R})_{i}{\bf e}_{i}$ may naively look like a scalar product, but actually
$({\bf w}_{R})_{i}$ are numerical coefficients multiplying the array ${\bf
e}_{i}$.
The dilatation along the flow, $D_{\parallel}$, is given by the projection
along ${\bf e}_{1}$:

$D_{\parallel}(t)={\bf D}(t)\cdot{\bf e}_{1}\simeq(\bm{\alpha}_{0}\cdot{\bf
w}_{R})\exp(\Lambda_{R}\,t)({\bf w}_{R})_{1}$ (22)

Since ${\bf w}_{R}$ is a unit vector, the dilatation perpendicular to the flow
is

$\displaystyle D_{\perp}(t)$ $\displaystyle=$ $\displaystyle\sqrt{({\bf
D}(t)\cdot{\bf e}_{2})^{2}+({\bf D}(t)\cdot{\bf e}_{3})^{2}}$ (23)
$\displaystyle\simeq$ $\displaystyle|\bm{\alpha}_{0}\cdot{\bf
w}_{R}|\exp(\Lambda_{R}\,t)\sqrt{({\bf w}_{R})_{2}^{2}+({\bf w}_{R})_{3}^{2}}$
$\displaystyle\equiv$ $\displaystyle|\bm{\alpha}_{0}\cdot{\bf
w}_{R}|\exp(\Lambda_{R}\,t)\sqrt{1-({\bf w}_{R})_{1}^{2}}$

We therefore conjecture the following test for chaotic motion: (i) one has
first of all to identify potentially unstable regions, i.e., where
$\Lambda_{R}>0$, and where therefore exponential dilatation of neighboring
trajectories occurs. (ii) Secondly, within these regions, one estimates $({\bf
w}_{R})_{1}$ as a proxy for the fraction of dilatation that takes place along
a direction parallel to the local flow. This dilatation does not contribute to
chaotic motion therefore, if it is dominant–i.e. $|({\bf w}_{R})_{1}|$ is
fairly close to unity–we may argue that dynamics is regular even in the
presence of positive eigenvalues. We will test this conjecture in the next
section.

##  III Analysis of some two-dimensional test cases

We consider here three study cases: (A) The Hènon-Heiles potential, and (B)
The anti-Henon-Heiles potential. They are classical cases, chosen in order to
study reference cases for chaotic and regular potential. (C) The third
potential was designed by ours, it is a chaotic potential as well.

If we write $p_{x}={\cal{K}}\cos\phi,\,p_{y}={\cal{K}}\sin\phi$, with $\cal K$
defined below Eq. (13), the eigenvalues $\Lambda_{R},\,\Lambda_{I}$ depend
upon three coordinates: $x,y,\phi$. We may envision two modalities of
operation: either we consider all the three coordinates $(x,y,\phi)$ as
independent parameters and scan the stability of the system within this three-
dimensional space. This is the same approach tackled by TBD but, in their
case, only the $x,y$ coordinates appeared within the analytical expression for
the rates $\mu$ in Eq. (8). We have to deal with a further degree of freedom,
which makes the scan more cumbersome and scarcely instructive. We prefer,
thus, the other way round. We produce a single trajectory by integrating
numerically the Hamilton equations (using a symplectic algorithm [27]). Then,
inspect the time series $(x(t_{i}),y(t_{i}),p_{x}(t_{i}),p_{y}(t_{i}))$ and
identify it be either regular or chaotic. Thirdly, we evaluate the relevant
quantities (eigenvalues and associated eigenvectors) over this time series.

###  III.1 The Hènon-Heiles system

The Henon-Heiles potential has the form

$V={1\over 2}(x^{2}+y^{2})+x^{2}y-{1\over 3}y^{3}$ (24)

This system is known to feature an order-chaos transition at energy $E\approx
0.083$, and to be totally chaotic when $E=1/6$ (above this value the
equipotential curves are no longer closed. We will consider here only bounded
systems) [6, 22, 26].

Figure 1 shows the results corresponding to a case with energy $E=0.12$: the
phase space is only partially chaotic, and the chosen trajectory appears
almost regular. The upper plot shows the $(x,y)$ motion of the numerical
trajectory: the black solid curve. The energetically allowed region is
surrounded by the boundary $V=E$. The regions unstable according to the TBD
criterion are filled with gray. Eigenvalues $\Lambda_{R}$ of the $\hat{\bf W}$
matrix are computed at different positions along the trajectory:
positive/negative eigenvalues $\Lambda_{R}$ are marked with red/green dots.
Remarkably, the trajectory does not cross any unstable region according to the
TBD criterion, whereas a fraction of it is flagged as unstable according to
Eq. (18). However, we have pointed out in the previous section that
instability does imply chaos only if there is appreciable dilatation along a
direction perpendicular to the trajectory. We attempt to illustrate this in
the lower plot, which features three time traces. The gray curve is the real-
valued eigenvalue $\Lambda_{R}$ from Eq. (18). The cyan curve is $({\bf
w}_{R})_{1}^{2}$. We have explained a few lines above that $({\bf w}_{R})_{1}$
makes sense as a marker of chaos only as long as $\Lambda_{R}>0$, therefore we
graphically emphasize this condition in the plot: the regions with
$\Lambda_{R}>0$ are those plotted with thicker and darker cyan tract. Finally,
the blue curve is the exponential rate of divergence according to the TBD
criterion, $\lambda_{TBD}$–which is here always negative. Throughout the whole
time trace the value of $|({\bf w}_{R})_{1}|$ is quite close to one, meaning
that most of the separation between trajectories occurs along the flow and
does not produce chaos. In Figure 2, finally, we show the corresponding time
traces of $\Lambda_{R},\Lambda_{I,r},\Lambda_{I,i}$. At least in this case, it
is apparent that only $\Lambda_{R}$ takes appreciably positive values.

Figure 1: An example of trajectory of the Hènon-Heiles system at $E=0.12$.
Upper plot: the $(x,y)$ projection of the orbit (black solid curve); the
energy boundary $V=E$ (gray curve). The gray-filled regions are those where
the TBD criterion predicts positive exponential rate of divergence. The green
and red dots superimposed to the particle trajectory flag the point where
$\Lambda_{R}<0$ (respectively $\Lambda_{R}>0$). Lower plot: gray curve is the
real-valued $\Lambda_{R}$ from Eq. (18); the cyan curve is $({\bf
w}_{R})_{1}^{2}$ plotted in correspondence of $\Lambda_{R}>0$ (darker and
thicker) or $\Lambda_{R}<0$ (lighter and thinner). The blue curve is the
eigenvalue computed according to the TBD criterion.  Figure 2: Time traces of
$\Lambda_{R},\Lambda_{I,r},\Lambda_{I,i}$ for the trajectory plotted in Figure
1.

For comparison, in Figure 3 we show another trajectory, at $E=0.15$. The
trajectory is now plainly irregular and, correspondingly, $|({\bf
w}_{R})_{1}|$ is noticeably lesser than unity.

Figure 3: Another example of trajectory of the Hènon-Heiles system, now at
$E=0.15$. Upper plot: the $(x,y)$ projection of the orbit (black solid curve);
the energy boundary $V=E$ (gray curve). The gray-filled regions are those
where the TBD criterion predicts positive exponential rate of divergence. The
green and red dots superimposed to the particle trajectory flag the point
where $\Lambda_{R}<0$ (respectively $\Lambda_{R}>0$). Lower plot: gray curve
is the real-valued $\Lambda_{R}$ from Eq. (18); the cyan curve is $({\bf
w}_{R})_{1}^{2}$ plotted in correspondence of $\Lambda_{R}>0$ (darker and
thicker) or $\Lambda_{R}<0$ (lighter and thinner). The blue curve is the
eigenvalue computed according to the TBD criterion.

###  III.2 The anti-Hènon-Heiles system

The potential differs from the previous case just by the sign of the cubic
term in $y$:

$V={1\over 2}(x^{2}+y^{2})+x^{2}y+{1\over 3}y^{3}$ (25)

It may be shown that this system is always integrable (see, e.g., [23]). The
TBD criterion, on the other hand, claims that this system is chaotic above
$E=1/24=0.04167$ (see [22]).

Figure 4: An example of trajectory of the anti-Hènon-Heiles system at
$E=0.062$. Upper plot: the $(x,y)$ projection of the orbit (black solid
curve); the energy boundary $V=E$ (gray curve). The gray-filled regions are
those where the TBD criterion predicts positive exponential rate of
divergence. The green and red dots superimposed to the particle trajectory
flag the point where $\Lambda_{R}<0$ (respectively $\Lambda_{R}>0$). Lower
plot: gray curve is the real-valued $\Lambda_{R}$ from Eq. (18); the cyan
curve is $({\bf w}_{R})_{1}^{2}$ plotted in correspondence of $\Lambda_{R}>0$
(darker and thicker) or $\Lambda_{R}<0$ (lighter and thinner). The blue curve
is the eigenvalue computed according to the TBD criterion.

Figure 4 shows an instance of the results. The trajectory overlaps an unstable
region according to TBD criterion, and positive and negative eigenvalues
computed according to the present method are spread along the whole
trajectory. Hence, inspection of the eigenvalues would naively to infer a
possible instability. Practically the whole divergence occurs parallel to the
flow, though: see the bottom plot, where $|({\bf w}_{R})_{1}|\approx 1$
everywhere, hence it does not impact upon the regularity of the motion.

###  III.3 A potential with mixed convex-concave energy boundary

As a third example, in order to enlarge the database, we choose the potential

$V={1\over 2}(x^{2}-y^{2})+a\exp(x^{2}+y^{2})$ (26)

This potential is interesting for a number of reasons. First of all, when
$a=0$ the system is integrable, yet it involves exponentially growing
solutions (and therefore also relative displacements) along the y-direction.
Secondly, when $a\neq 0$ its isopotential contour levels become closed but,
unlike most of the potentials commonly encountered, they are not everywhere
convex (bulging outwards), rather alternating convex and concave intervals
exist. The role of the energy boundary upon the regular or chaotic kind of
motion of a system has often been taken into consideration since it is the
region where the strongest acceleration does take place. Billiards are a class
of dynamical systems characterized by step-like potentials. It is well known
that the shape of the boundary determines their chaotic or regular nature, and
that billiards having a concave fraction of the boundary are chaotic [28].
Furthermore, to some extent the properties of the billiards are inherited by
smooth dynamical systems [29]. Rod, Pecelli and Churchill [30] showed that
periodic orbits are hyperbolically unstable in smooth potentials whose
boundary has concave components alongside some (relatively mild) symmetry
requirements. A similar result was obtained more recently by different methods
by Li and Zhang [31].
In the following, we show results for $a=2\times 10^{-4},E=0.5$. The Poincaré
plot (5) ensures that, for these parameters, the motion is actually chaotic.
Figure 6 is the counterpart of figures (1,3,4). This time, the oscillations
between the parallel and perpendicular directions are wild, qualitatively we
may claim that this system appears more chaotic than the Hènon-Heiles one.

Figure 5: Poincarè plot for a trajectory with the potential (26), obtained
taking the intersections with the plane $y=0$ and ${\dot{y}}<0$. Figure 6:  An
example of trajectory for the system defined by the potential Eq. (26), with
$a=2\times 10^{-4},E=0.5$. Upper plot: the $(x,y)$ projection of the orbit
(black solid curve); the energy boundary $V=E$ (gray curve). The gray-filled
regions are those where the TBD criterion predicts positive exponential rate
of divergence. The green and red dots superimposed to the particle trajectory
flag the point where $\Lambda_{R}<0$ (respectively $\Lambda_{R}>0$) . Lower
plot: gray curve is the real-valued $\Lambda_{R}$ from Eq. (18); the cyan
curve is $({\bf w}_{R})_{1}^{2}$ plotted in correspondence of $\Lambda_{R}>0$
(darker and thicker) or $\Lambda_{R}<0$ (lighter and thinner). The blue curve
is the eigenvalue computed according to the TBD criterion.

##  IV Conclusions

The shortcoming of the TBD criterion is that it treats on the same par couples
of trajectories that are informative about the dynamics of the system, and
trajectories that are not, either because do not lie both on the same energy
surface, or because are both aligned along the common flow. In this work we
have shown how–at least in the two-dimensional case–one can implement an
amended criterion: chaotic systems are those which feature a positive
exponential rate of dilatation along the constant-energy hypersurface when the
dilatation itself is appreciably aligned along directions perpendicular to the
local flow. We did not establish a threshold quantifying what an appreciable
alignment means, however it appears intuitive in the cases shown. The proposed
criterion turns out to provide sensible interpretations of the numerical
results in all the test cases considered. Unlike the original unreliable TBD
criterion–and like the other criteria listed in the Introduction–it relies
upon the inspection of pre-computed trajectories, hence cannot be employed as
a shortcut to the actual solving of Hamilton’s equations.

## Acknowledgments

Barbara Momo read carefully a draft of this work, making several useful
suggestions. F.S. acknowledges fruitful discussions with D.F. Escande and I.
Predebon. L.S. is partially supported by the European Quantum Flagship Project
”PASQuanS 2”, by the European Union-NextGenerationEU within the National
Center for HPC, Big Data and Quantum Computing [Project No. CN00000013, CN1
Spoke 10: “Quantum Computing”], by the BIRD Project ”Ultracold atoms in curved
geometries” of the University of Padova, and by “Iniziativa Specifica Quantum”
of INFN.

## References

  * [1] C. Skokos, J. Phys. A 34, 10029 (2001)
  * [2] C.L.Vozikis, H. Varvoglis, K. Tsiganis, A&A 359, 386 (2000)
  * [3] C. Froeschlé, E. Lega and R. Gonzi R., Celest. Mech. Dyn. Astron. 67, 41 (1997)
  * [4] J.P. Codaccioni, F. Doveil, D.F. Escande, Phys. Rev. Lett. 49, 1879 (1982)
  * [5] A. Okubo, Deep Sea Res. 17, 445 (1970)
  * [6] M. Toda, Phys. Lett. 48A, 335 (1974)
  * [7] P. Brumer and J. Duff, J. Chem. Phys. 65, 3566 (1976)
  * [8] C.P. Enz, M.O. Hongler and C.V. Quach Thi, Helv. Phys. Acta 48, 787 (1975)
  * [9] H.N. Núnez-Yépez, A.L. Salas-Brito, C.A. Vargas, L. Vicente, Phys. Lett. A 145, 101 (1990)
  * [10] A. Oloumi, D. Teychenné, Phys. Rev. E 60, R6279 (1999)
  * [11] T. Watabe, Phys. Lett. B 343, 254 (1995)
  * [12] V.I. Kuvnishov and A.V. Kuzmin, Phs. Lett. A 296 , 82 (2002)
  * [13] V.I. Kuvnishov and V.A. Piatrou, Phys. Part. Nuclei Lett. 5, 72 (2008)
  * [14] Li Junqing, Zhu Jieding, Gu Jinnan, Phys. Rev. B 52, 6458 (1995)
  * [15] V.R. Manfredi and L. Salasnich, Mod. Phys. Lett. A 12, 1851 (1997)
  * [16] P.S. Letelier, J. Ramos-Caro, F. López-Suspes, Phys. Lett. A 375, 3655 (2011)
  * [17] A. Ghosh et al, Chaos 24, 013117 (2014)
  * [18] Yu-L. Chang, Lie-Y. Oey, Ocean Dyn. 64, 259 (2014)
  * [19] F.J. Beron-Vera, A.Hadjighasem, Q. Xia, M.J. Olascoaga, and G. Haller, PNAS 116, 18251 (2019)
  * [20] B.K. Shivamoggi, G.J.F. van Heijst and L.P.J. Kamp, Fluid Dynamics Res 54, 015505 (2022)
  * [21] G. Benettin, R. Brambilla, L. Galgani, Physica 87A, 381 (1977)
  * [22] A.K. Pattanayak and W.C. Schieve, Z. Naturforsch. 52a, 34 (1997)
  * [23] R.G. Smirnov, Appl. Math. Lett. 11, 71 (1998)
  * [24] C. Cerjan, W.P. Reinhardt, J. Chem. Phys. 71, 1819 (1979)
  * [25] R. Kosloff and S.A. Rice, J. Chem. Phys. 74, 1947 (1981)
  * [26] E.E. Zotos, Nonlinear Dyn. 79, 1665 (2015)
  * [27] R.I. McLachlan and P. Atela, Nonlinearity 5, 541 (1992)
  * [28] Ya.G. Sinai, Russian Math. Surveys 25, 137 (1970)
  * [29] A. Rapoport, V. Rom-Kedar, D. Turaev, Commun. Math. Phys. 272, 567 (2007)
  * [30] D.L. Rod, G. Pecelli, R.C. Churchill, Jour. Diff. Eq. 24, 329 (1977)
  * [31] J. Li, S. Zhang, Phys. Lett. A 375, 1710 (2011)