# Warming up recurrent neural networks to maximise reachable multistability
greatly improves learning

Gaspard Lambrechts
<EMAIL_ADDRESS>&Florent De Geeter
<EMAIL_ADDRESS>&Nicolas Vecoven
<EMAIL_ADDRESS>&Damien Ernst
<EMAIL_ADDRESS>&Guillaume Drion
<EMAIL_ADDRESS>
###### Abstract

Training recurrent neural networks is known to be difficult when time
dependencies become long. Consequently, training standard gated cells such as
the gated recurrent unit (GRU) and the long short-term memory (LSTM) on
benchmarks where long-term memory is required remains an arduous task. In this
work, we show that although most classical networks have only one stable
equilibrium at initialisation, learning on tasks that require long-term memory
only occurs once the number of network stable equilibria increases; a property
known as multistability. Multistability is often not easily attained by
initially monostable networks, making learning of long-term dependencies
difficult. This insight leads to the design of a novel, general way to
initialise any recurrent network connectivity through a procedure called
“warmup” to improve its capability to learn arbitrarily long time
dependencies. This initialisation procedure is designed to maximise network
reachable multistability, i.e., the number of attractors within the network
that can be reached through relevant input trajectories. Warming up is
performed before training, using stochastic gradient descent on a specifically
designed loss. We show on information restitution, sequence classification,
and reinforcement learning benchmarks that warming up greatly improves
recurrent neural network performance for multiple recurrent cell types, but
sometimes impedes precision. We therefore introduce a parallel recurrent
network structure with a partial warmup that is shown to greatly improve
learning of long-term dependencies in sequences while maintaining high levels
of precision. This approach provides a general framework for improving
learning abilities of any recurrent cell type when long-term memory is
required.

_Keywords_ recurrent neural network $\cdot$ multistability $\cdot$
initialisation procedure $\cdot$ long-term memory $\cdot$ warmup

##  1 Introduction

Despite their performances and widespread use, recurrent neural networks
(RNNs) have been known to be blackbox models with extremely complex internal
dynamics. A growing body of work has focused on understanding the internal
dynamics of trained RNNs [Sussillo and Barak, 2013, Ceni et al., 2020,
Maheswaranathan et al., 2019], providing invaluable insights into the RNN
prediction process. In particular, recent work has highlighted the important
role played by fixed points in RNNs, and it has argued that locating such
fixed points efficiently could provide insights into RNN dynamics and input-
output properties [Katz and Reggia, 2017]. Here, we build upon this line of
work to highlight how maximizing the number of reachable fixed points before
learning improve RNN learning, in particular in the presence of arbitrarily
long-term dependencies.

Training RNNs has been shown to be difficult when time dependencies become too
long [Pascanu et al., 2013]. This limitation is often due to vanishing or
exploding gradients. Different methods have been proposed to help with this
issue, such as introducing gating mechanisms [Cho et al., 2014, Hochreiter and
Schmidhuber, 1997], clipping gradients [Pascanu et al., 2013] and also
maintaining orthogonality in recurrent weight matrices [Jing et al., 2019].
Other authors have also proposed new ways of initialising RNN parameters
[Pascanu et al., 2013, Van Der Westhuizen and Lasenby, 2018, Marichal et al.,
2009], endowing them with better training properties. Likewise, novel
recurrent cell dynamics, such as neuromodulated bistable recurrent cells
(nBRCs) [Vecoven et al., 2021], have been introduced to help tackle long-term
memory benchmarks. nBRCs were specifically designed to maximise reachability
of cellular bistability, providing a way to create never-fading memory at the
cellular level. These results highlighted how dynamics of untrained RNNs,
i.e., at initialisation, can strongly impact learning performance of RNNs.

In this work, we extend this approach at the network level by maximizing
reachable multistability of any recurrent cell type prior to learning. To this
end, we propose a novel RNN pre-training procedure called “warmup” that is
designed to maximise the number of RNN attractors that can be reached by a set
of training samples. First, we introduce a fast-to-compute measure called
variability amongst attractors (VAA) that counts the number of different
reachable attractors within a network. We show that loss decrease during
learning in long-term memory benchmarks is highly correlated with an increase
in VAA, highlighting both the relevance of the measure and the importance of
multistability for efficient learning. Second, we use stochastic gradient
ascent on a diffentiable proxy of the VAA before training as a way of
maximizing the number of reachable attractors within the network. We show that
this technique strongly improves performance on long-term memory benchmarks,
at the cost of precision, the latter relying on the richness of network
transient dynamics. Finally, we propose a parallel recurrent network structure
with a partial warmup that enables one to combine long-term memory through
multistability with precision through rich transient dynamics. For multiple
RNNs such as gated recurrent units (GRUs) [Cho et al., 2014] and long short-
term memories (LSTMs) [Hochreiter and Schmidhuber, 1997], we show that this
method indeed retains the benefits of the warmup while increasing predictive
performances.

In Section 2, RNNs are introduced as dynamical systems and the concept of
multistability is introduced for those systems. In Section 3, the supervised
learning and reinforcement learning benchmarks are introduced. In Section 4,
the VAA is introduced as a measure of the multistability given initial states
and its estimation procedure is detailed for states resulting from a set of
input sequences. The correlation between multistability and learning is shown
empirically on the benchmarks. In Section 5, a differentiable proxy of the VAA
is introduced and the warmup procedure that fosters multistability at
initialisation is given. The benefits of warmup are shown empirically on long-
term information restitution benchmarks. In addition, the double-layer
architecture is introduced to achieve a better performance both on information
restitution benchmarks and sequence classification benchmarks. Finally,
Section 6 concludes and proposes several future works.

##  2 Background

In this section, RNNs are introduced as dynamical systems. The fixed points of
these systems are defined, and the notions of attractors, reachable attractors
and multistability are introduced.

###  2.1 Recurrent neural networks

RNNs are parametric function approximators that are often used to tackle
problems with temporal structure. Indeed, RNNs process the inputs
sequentially, exhibiting memory through hidden states that are outputted after
each timestep, and processed at the next timestep along with the following
input. These connections allow RNNs to memorise relevant information that
should be captured over multiple timesteps. More formally, an RNN architecture
is defined by its update function $f$, its output function $g$ and its
initialisation function $h$ that are parameterised by a parameter vector
$\theta\in\mathbb{R}^{d}$. Let
$\mathbf{u}_{1:T}=[\mathbf{u}_{1},\dots,\mathbf{u}_{T}]$, with
$T\in\mathbb{N}$ and $\mathbf{u}_{t}\in\mathbb{R}^{n}$, an input sequence.
RNNs maintain an internal memory state $\mathbf{x}_{t}$ through an update rule
$\mathbf{x}_{t}=f(\mathbf{x}_{t-1},\mathbf{u}_{t};\theta)$ and output a value
$\mathbf{o}_{t}=g(\mathbf{x}_{t};\theta)$, where the initial hidden state
$\mathbf{x}_{0}=h(\theta)$ is often chosen to be zero. We note that often, the
output of the RNN is simply its hidden state $\mathbf{x}_{t}$, i.e. $g$ is the
identity function. RNNs can be composed of only one recurrent layer, or they
can be built with multiple layers that are linked sequentially through
$\mathbf{u}^{i}_{t}=\mathbf{o}^{i-1}_{t}$ with
$\mathbf{u}^{0}_{t}=\mathbf{u}_{t}$, where $\mathbf{o}^{i}_{t}$ denotes the
output of layer $i$ and $\mathbf{u}^{i}_{t}$ its input. Backpropagation
through time is used to train these networks where gradients are computed
through the complete sequence via the hidden states [Werbos, 1990]. In
practice, it is difficult to do this on longer sequences, due to vanishing and
exploding gradients. For an in-depth study on the subject and some proposed
solutions, we refer the reader to [Pascanu et al., 2013]. The following
recurrent architectures are considered in the experiments: long short-term
memory (LSTM) [Hochreiter and Schmidhuber, 1997], gated recurrent unit (GRU)
[Chung et al., 2014], bistable recurrent cell (BRC) and recurrently
neuromodulated BRC (nBRC) [Vecoven et al., 2021], and minimal gated unit (MGU)
[Zhou et al., 2016]. In addition, we consider the chrono-initialised LSTMs
[Tallec and Ollivier, 2018].

###  2.2 Fixed points in recurrent neural networks

In dynamical systems, fixed points are defined as points in the phase space
that map to themselves through the update function, for a given input
$\mathbf{u}$. For a system $f$, we say that a state $\mathbf{x}^{*}$ is a
fixed point in $\mathbf{u}$ if and only if

$\mathbf{x}^{*}=f(\mathbf{x}^{*},\mathbf{u}).$ (1)

Fixed points can either be fully attractive (attractors), fully repulsive
(repellors), or combine attractive and repulsive manifolds (saddle points).
The set of starting states for which the system converges to the fixed point
$\mathbf{x}^{*}$ is called basin of attraction of $\mathbf{x}^{*}$ in
$\mathbf{u}$ and is written as

$\displaystyle\mathcal{B}^{\mathbf{u}}_{\mathbf{x}^{*}}$
$\displaystyle=\left\\{\mathbf{x}\;\middle|\;\lim_{n\to\infty}f^{n}(\mathbf{x},\mathbf{u})=\mathbf{x}^{*}\right\\}$
(2) $\displaystyle\text{ with }f^{n}(\mathbf{x},\mathbf{u})$
$\displaystyle=\begin{cases}f\left(f^{n-1}(\mathbf{x},\mathbf{u}),\mathbf{u}\right)\text{
if $n>1$},\\\ f(\mathbf{x},\mathbf{u})\text{ if $n=1$}.\end{cases}$ (3)

If the limit is not defined for some point $\mathbf{x}$, then this point does
not belong to any basin of attraction in $\mathbf{u}$. No attractor will ever
be reached by starting from this point. Mathematically, $\mathbf{x}^{*}$ is an
attractor in $\mathbf{u}$ if its basin of attraction in $\mathbf{u}$,
$\mathcal{B}^{\mathbf{u}}_{\mathbf{x}^{*}}$, is not a singleton.

From this definition, one can introduce the notion of reachable attractors. In
particular, we say that an attractor $\mathbf{x}^{*}$ in $\mathbf{u}$ is
reachable for some initial state $\mathbf{x}_{0}$ and input sequence
$\mathbf{u}_{1:T}$ if and only if
$\mathbf{x}_{T}\in\mathcal{B}^{\mathbf{u}}_{\mathbf{x}^{*}}$, with
$\mathbf{x}_{t}=f(\mathbf{x}_{t-1},\mathbf{u}_{t}),\;t=1,\dots,T$. Given an
initial state and a batch of input sequences, a system that has a single
reachable attractor in $\mathbf{u}$ for all input sequences is said to be
monostable in $\mathbf{u}$ for the batch of input sequences, whereas a system
that has multiple reachable attractors in $\mathbf{u}$ is said to be
multistable in $\mathbf{u}$ for the batch of input sequences.

Due to their temporal nature and update rules, RNNs can be seen as discrete-
time non-linear dynamical systems. Formally, given a parameter vector
$\theta$, the system $f$ is given by the update function of the RNN, such that
$f(\mathbf{x},\mathbf{u})=f(\mathbf{x},\mathbf{u};\,\theta)$. This viewpoint
has already been used to understand the difficulties for RNNs to capture
longer-time dependencies [Bengio et al., 1993, Doya, 1993]. Recently, it was
shown that fixed points are important when it comes to understanding the
prediction process of trained RNNs [Sussillo and Barak, 2013]. Since
attractors correspond to network steady states, they are thought to be the
allowing factor for RNNs to retain information over a long period of time
[Pascanu et al., 2013, Sussillo and Barak, 2013, Maheswaranathan et al.,
2019].

##  3 Benchmarks

In this section, the different benchmarks are introduced. First, the
supervised learning tasks are introduced, including long-term information
restitution benchmarks in Subsection 3.1 and sequence classification
benchmarks in Subsection 3.2. In Subsection 3.3, a reinforcement learning
benchmark with partially observable environment is introduced. This
environment contains long-term dependencies.

###  3.1 Long-term information restitution benchmarks

The benchmarks introduced in this subsection contain long-term dependencies,
and therefore require networks able to remember relevant information for a
long period. These will be referred to as long-term information restitution
benchmarks, or in shorter form, information restitution benchmarks.

#### Copy first input benchmark.

In this benchmark, the network is presented with a one-dimensional sequence of
$T$ timesteps $\mathbf{u}_{1:T}$, where
$\mathbf{u}_{t}\sim\mathcal{N}(0,1),\;t=1,\dots,T$, and is tasked at
approximating the target $\mathbf{y}_{T}=\mathbf{u}_{1}$. This benchmark thus
consists of memorizing the initial input for $T$ timesteps. It allows one to
measure the ability of recurrent architectures to bridge long-term
dependencies when the length $T$ is large. Given the output $\mathbf{o}_{T}$
of the network, we seek to minimise the squared error
$\mathcal{L}(\mathbf{o}_{T},\mathbf{y}_{T})=(\mathbf{o}_{T}-\mathbf{y}_{T})^{2}$.

#### Denoising benchmark.

In this benchmark, the network is presented with a two-dimensional sequence of
$T$ timesteps. The first dimension is a noised input stream
$\mathbf{u}_{1:T}^{1}$, where
$\mathbf{u}_{t}^{1}\sim\mathcal{N}(0,1),\;t=1,\dots,T$. Five timesteps of this
stream should be remembered and outputted one by one by the network at
timesteps $\\{T-4,\dots,T\\}$. These five timesteps
$\mathcal{S}=\\{t_{1},t_{2},t_{3},t_{4},t_{5}\\}$, with
$t_{1}<t_{2}<t_{3}<t_{4}<t_{5}$, are sampled without replacement in
$\\{1,\dots,T-N\\}$ with $N\geq 5$. $N$ is a hyperparameter that allows one to
tune how long the network should be able to retain the information at a
minimum. The five timesteps are communicated to the network through the second
dimension of the input $\mathbf{u}_{1:T}^{2}$, where $\mathbf{u}_{t}^{2}=1$ if
$t\in\mathcal{S}$, and $\mathbf{u}_{t}^{2}=0$ otherwise, for $t=1,\dots,T$.
The target is thus given by
$\mathbf{y}_{T-4:T}=[\mathbf{u}_{t_{1}},\mathbf{u}_{t_{2}},\mathbf{u}_{t_{3}},\mathbf{u}_{t_{4}},\mathbf{u}_{t_{5}}]$.
Given the output sequence $\mathbf{o}_{T-4:T}$ of the network, we seek to
minimise the mean squared error
$\mathcal{L}(\mathbf{o}_{T-4:T},\mathbf{y}_{T-4:T})=\sum_{t=T-4}^{T}(\mathbf{o}_{t}-\mathbf{y}_{t})^{2}$.

###  3.2 Sequence classification benchmarks

The benchmarks introduced in this subsection are sequence classification
problems and therefore require networks able to use the information received
in the sequence in order to infer the class. These will be referred to as
_sequence classification benchmarks_.

#### Permuted sequential MNIST.

In this benchmark, the network is presented with the MNIST images, where
pixels are presented to the network one by one as a sequence of length
$T=28\times 28=784$. It differs from the regular sequential MNIST in that
pixels are shuffled in a random order. Note that all images are shuffled
according to the same random order. The network is tasked at outputting a
probability for each possible digit that could be represented in the initial
image. This benchmark is known to be a more complex challenge than the regular
one. Given the output $\mathbf{o}_{T}\in\mathbb{R}^{10}$ of the network and
the true digit index $\mathbf{y}_{T}\in\\{1,\dots,10\\}$, we seek to minimise
the negative log likelihood loss
$\mathcal{L}(\mathbf{o}_{T},\mathbf{y}_{T})=-\log(\mathbf{o}_{T}^{\mathbf{y}_{T}})$.

#### Permuted line-sequential MNIST.

This benchmark is the same as the permuted sequential MNIST benchmark, except
that the pixels are fed $28$ by $28$, which corresponds to one line of the
permuted image. The input dimension is thus $28$ instead of one. $N$ black
lines are added at the end of the sequence such that the total length of the
sequence is $T=28+N$. This has the effect of a forgetting period, such that
any relevant information for predicting the class probabilities will be
farther from the prediction timestep $T$.

###  3.3 Reinforcement learning benchmark

In reinforcement learning, the function approximators also process sequences
as input when considering partially observable Markov decision processes
(POMDPs). Indeed, in such environments, the optimal policies, as well as the
value functions, are functions of the complete sequence of observations and
past actions, called the history. In this work, we focus on the approximation
of the history-action value function, or $\mathcal{Q}$-function, in order to
derive a near-optimal policy in the considered POMDP. The deep recurrent
Q-network (DRQN) algorithm is used to approximate this $\mathcal{Q}$-function
with an RNN. From this approximation, we derive the fully greedy policy by
taking the action that maximises the $\mathcal{Q}$-function for any given
history. See A for the formal definition of POMDPs and their
$\mathcal{Q}$-functions, and see B for the detailed DRQN algorithm.

Figure 1: T-Maze layout example, with the initial position of the agent in
black, the treasure in green and the cell to avoid in red.

The partially observable environment that is considered is the T-Maze
environment [Bakker, 2001]. The T-Maze is a POMDP where the agent is tasked
with finding the treasure in a T-shaped maze (see Figure 1). The state is
given by the position of the agent in the maze and the maze layout that
indicates whether the goal lies up or down after the crossroads. The initial
state determines the maze layout, and it never changes afterwards. The initial
observation made by the agent indicates the layout. Navigating in the maze
provides zero reward, except when bouncing onto a wall, in which case a reward
of $-0.1$ is received. While traveling along the maze, the agent only receives
the information that it has not yet reached the junction. Once the junction
reached, the agent is notified: it must now choose a direction depending on
the past information it remembers. Finding the treasure provides a reward of
$4$. Passed the crossroads, the states are always terminal. The optimal policy
thus consists of going through the maze, while remembering the initial
observation in order to take the correct direction at the crossroads. This
POMDP is parameterised by the corridor length $L\in\mathbb{N}$ that determines
the number of timesteps for which the agent should remember the initial
observation. The discount factor is $\gamma=0.98$. This POMDP is formally
defined in C.

##  4 Correlating multistability and learning

This section aims at showing the correlation that exists between
multistability properties of RNNs and their ability to learn long-term
dependencies. To this end, in Subsection 4.1 we first introduce the VAA, a
measure of the number of basins of attraction that are spanned by a set of
states. In Subsection 4.2, we show how to estimate the multistability of an
RNN using VAA by computing the number of reachable attractors for a batch of
input sequences. We then carry out a number of experiments in Subsection 4.3
to show the correlation between multistability and learning with different
types of RNN on the benchmarks previously introduced.

###  4.1 Variability amongst attractors

One way to measure the multistability of a system is to count the number of
different attractors that can be reached starting from different initial
states. Our measure is named Variability Amongst Attractors (VAA). The VAA of
a system $f$ for a set of initial states $\mathcal{X}$ and an input
$\mathbf{u}$ is defined as

$\text{VAA}(f,\mathcal{X},\mathbf{u})=\frac{1}{\left\lvert\mathcal{X}\right\rvert}\sum^{\left\lvert\mathcal{X}\right\rvert}_{i=1}\frac{1}{\sum\limits_{j=1}^{\left\lvert\mathcal{X}\right\rvert}\delta\left(\left\lVert\lim\limits_{n\to\infty}f^{n}(\mathbf{x}_{i},\mathbf{u})-\lim\limits_{n\to\infty}f^{n}(\mathbf{x}_{j},\mathbf{u})\right\lVert=0\right)}$
(4)

where $\delta(x)$ is the Kronecker delta function that returns $1$ when
condition $x$ is met, and $0$ otherwise. It is assumed that the system $f$ is
such that limits exist.

The denominator of Equation 4 gives the number of states in $\mathcal{X}$ that
converge towards the same attractor as $\mathbf{x}_{i}$. The sum of this
fraction over all states that converge towards a given attractor is thus equal
to one, such that the sum of this fraction over all states gives the number of
different attractors. $\text{VAA}(f,\mathcal{X},\mathbf{u})$ is thus equal to
the number of different attractors in $\mathbf{u}$ reached from the initial
states contained in $\mathcal{X}$ divided by the number of initial states
$\left\lvert\mathcal{X}\right\rvert$. Its maximal value is thus $1$, when all
reached attractors are different, and its minimal value is
$\frac{1}{\left\lvert\mathcal{X}\right\rvert}$, when all the states have
converged towards the same attractor (i.e., the system is monostable).

###  4.2 Using the VAA to measure the multistability of an RNN

RNNs can exhibit a long-lasting memory through multistability in their hidden
states [Vecoven et al., 2021]. Indeed, having multiple basins of attraction
reachable from different input sequences allows one to encode information
about these sequences over the long term. In this section, we propose
estimating the proportion of different basins of attraction that are reached
in an RNN for a batch of input sequences, using the VAA measure.

In practice, since it is impossible to evaluate the limits to infinity in the
VAA, we fix a finite number of timesteps $M$ for state convergence, called the
stabilization period. As a consequence, the system may not have completely
converged towards the attractor after this period. We thus define a tolerance
$\varepsilon$ below which two final states are considered to correspond to the
same attractor. In all experiments, $\varepsilon$ is chosen to
$1\text{\times}{10}^{-4}$. This truncated VAA is written as

$\text{VAA}_{M,\varepsilon}(f,\mathcal{X},\mathbf{u})=\frac{1}{\left\lvert\mathcal{X}\right\rvert}\sum^{\left\lvert\mathcal{X}\right\rvert}_{i=1}\dfrac{1}{\overset{\left\lvert\mathcal{X}\right\rvert}{\underset{j=1}{\sum}}\,\delta\left(\left\lVert
f^{M}(\mathbf{x}_{i},\mathbf{u})-f^{M}(\mathbf{x}_{j},\mathbf{u})\right\lVert\leq\varepsilon\right)}.$
(5)

In order to estimate the multistability of a network by its proportion of
reachable attractors for a given dataset, we propose estimating the mean value
of the VAA over several batches and input perturbations. To compute each value
of the VAA, a batch of sequences is sampled from the dataset, then a hidden
state is sampled at a random timestep for each sequence. This set of hidden
states constitutes the set of initial states used by Equation 5. Finally, we
sample a random perturbation
$\mathbf{u}\sim\mathcal{N}(\mathbf{0},\mathbf{1})$ that follows a standard
normal distribution. The proportion of attractors is then estimated with
$\text{VAA}_{M,\varepsilon}$. The resulting procedure is given in Algorithm 1.
In the experiments, the estimation is performed frequently throughout the
training procedure, but with a single batch for efficiency.

1

2

Parameters : $I\in\mathbb{N}$ the number of iterations to compute the mean of
the VAA.

$M\in\mathbb{N}$ the stabilization period.

$\varepsilon\in\mathbb{R}^{+}$ tolerance when considering state similarity.

$\theta\in\mathbb{R}^{d_{\theta}}$ the parameters of the network.

3

Data:
$\mathcal{D}=\\{\mathbf{u}_{1:T_{1}}^{1},\dots,\mathbf{u}_{1:T_{n}}^{N}\\}$ a
set of $N$ input sequences.

4

5Let $f\leftarrow f(\cdot,\cdot;\;\theta)$ the dynamical system.

6 Initialise mean value $\overline{\text{VAA}{}}\leftarrow 0$.

7 for  _$i=1,\dots,I$_ do

8 Sample a batch of input sequences $\mathcal{B}\sim\mathcal{D}$.

9

10 Sample a random hidden state in each input sequence
$\mathcal{X}\leftarrow\textnormal{{RandomHiddenStates}}(\mathcal{B},\theta)$.

11

12 Sample $\mathbf{u}\sim\mathcal{N}(\mathbf{0},\mathbf{1})$

13

14
$\overline{\text{VAA}{}}\leftarrow\overline{\text{VAA}{}}+\frac{1}{I}\,\text{VAA}{}_{M,\varepsilon}(f,\mathcal{X},\mathbf{u})$

15

16return _$\overline{\text{VAA}{}}$_

Algorithm 1 Estimating the proportion of reachable attractors of an RNN for a
set of input sequences

1

2

Parameters : $\theta\in\mathbb{R}^{d_{\theta}}$ the parameters of the network.

3

Data:
$\mathcal{B}=\\{\mathbf{u}_{1:T_{1}}^{1},\dots,\mathbf{u}_{1:T_{n}}^{n}\\}$ a
batch of $n$ input sequence sampled in the training set.

4

5$\mathcal{X}\leftarrow\\{\\}$

6

7foreach  _$\mathbf{u}_{1:T_{i}}^{i}\in\mathcal{B}$ _ do

8

9 Sample a timestep $t\sim\mathcal{U}(\left\\{1,\dots,T_{i}\right\\})$

10

11 Set $\mathbf{x}^{i}_{0}=h(\theta)$ where $h$ is the RNN’s initialisation
function

12

13 for  _$k=1,\dots,t$_ do

14 Set
$\mathbf{x}^{i}_{k}=f(\mathbf{x}^{i}_{k-1},\mathbf{u}_{k}^{i};\,\theta)$ where
$f$ is the RNN’s update function.

15

16 Update
$\mathcal{X}\leftarrow\mathcal{X}\bigcup\left\\{\mathbf{x}^{i}_{t}\right\\}$

17

18return _$\mathcal{X}$_

Algorithm 2 RandomHiddenStates

###  4.3 Experiments

In this subsection, we observe how the multistability of RNNs evolves when
they are trained on the long-term information restitution and reinforcement
learning benchmarks introduced in Section 3. The multistability of these
networks is estimated throughout the training procedure using Algorithm 1.

#### Copy first input benchmark

Figure 2 shows the performance of the different cells on this benchmark for
different sequence lengths $T\in\\{50,300,600\\}$. The best-performing cells
are BRC and nBRC, whose performances are not affected by the length of the
sequences. In comparison, the classical cells, MGU, LSTM and GRU, struggle to
decrease their losses. Furthermore, the longer the sequences are, the worse
their performances are. The last cell, the chrono-initialised LSTM, competes
with BRC and nBRC when its hyperparameter $T_{max}=600$ matches the length of
the sequences. When this is not the case, it still achieves good results, but
performs worse than BRC and nBRC. Figure 3 illustrates the correlation between
the VAA and the validation loss for the LSTM and CHRONO cells. The LSTM cell,
whose VAA never increases, does not manage to decrease its loss. On the other
hand, the CHRONO cell sees its loss decreasing while its VAA increases. This
figure also shows that the chrono initialisation promotes the learning of
long-term dependencies through multistability.

Figure 2: Mean square error ($\pm$ standard deviation) of different
architectures on the test set for the copy first input benchmark. Results are
shown after $50$ epochs and for different values of $T$.

Figure 3: Evolution of the validation loss (left) and of the VAA (right) of
LSTM networks, with and without chrono initialisation, trained on the copy
first input benchmark. The stabilization period for the VAA is chosen to
$M=$2000$$. The average and standard deviations over three runs are reported.

#### Denoising benchmark

Figure 4 shows the evolution of the VAA and the validation loss of multiple
GRU and MGU cells during the training on this benchmark. It is observed that
both types of network start to learn, i.e., their loss starts to decrease when
their VAA starts to increase. Figure 5 shows the validation loss and the VAA
of three different trainings of the CHRONO cell. Only the network that manages
to increase its VAA and maintain it at a high level achieves a good
performance. This proves once more the correlation between the VAA and the
learning on long-term information restitution benchmarks.

Figure 4: Evolution of the validation loss (left) and of the VAA (right) of
GRU and MGU networks trained on the denoising benchmark. The stabilization
period for the VAA is chosen to $M=$2000$$. The average and standard
deviations over three runs are reported. Learning, represented by a decreasing
loss, only starts when the network becomes multistable (VAA greater than
$\frac{1}{\left\lvert S\right\rvert}$).

Figure 5: Evolution of the validation loss (left) and of the VAA (right) of
multiple LSTM chrono networks trained on the denoising benchmark. The
stabilization period for the VAA is chosen to $M=$2000$$. Only the network
that managed to increase its VAA and to maintain it at a stable high value
learned to perform well.

#### T-Maze benchmark

In this reinforcement learning setting, a policy is derived from the
approximation of the $\mathcal{Q}$-function. On the left in Figure 6, we can
see the mean non-discounted cumulative reward obtained by the policies derived
from GRU cells approximating the $\mathcal{Q}$-function. On the right in
Figure 6, we can see the VAA of these cells estimated with Algorithm 1 using
the histories of the replay buffer as input sequences. Those value are clearly
correlated. Indeed, the better the agent plays, the higher its VAA is.

Figure 6: Evolution of the mean cumulative reward (left) obtained by GRU
agents and their VAA (right) during their training on a T-Maze of length 200.
The stabilization period for the VAA is chosen to $M=$10\,000$$. The average
and standard deviations over three runs are reported.

Figure 7: Evolution of the mean cumulative reward obtained by agents during
their training (left) and mean number of episodes required to reach the
optimal policy (right) on T-Mazes of length 20 (top) and 200 (bottom).

##  5 Fostering multistability at initialisation

In Subsection 5.1, we describe the warmup initialisation procedure that allows
one to maximise the multistability of a network for a dataset of input
sequences. Then, in Subsection 5.2, we compare classic cells to warmed-up
cells on information restitution, sequence classification, and RL benchmarks
and show the benefits of the warmup in tasks with long-term dependencies.
However, we also show that the warmup procedure does not improve the results
in the sequence classification tasks. Finally, in Subsection 5.3, we introduce
the double-layer architecture, that has both multistable and transient
dynamics. We show that this architecture reaches a better performance both on
information restitution and sequence classification benchmarks.

###  5.1 Warming up RNNs

The previous observations, that show a correlation between the multistability
of a network and its ability to learn long-term dependencies, suggest that
fostering multistability could ease learning in the case of tasks with long-
term dependencies. For instance, Figure 7 shows that multistable cells, such
as the BRC and nBRC cells, learn quickly on the T-Maze environment and obtain
the optimal policy faster than classical cells. In order to promote the
multistability of a network, we propose maximizing the number of reachable
attractors for random batches of input sequences coming from the dataset.

In practice, we propose using stochastic gradient descend (SGD) to maximise
the proportion of reachable attractors for batches of input sequences. We note
however that SGD cannot be used directly on the estimation of the proportion
of reachable attractors detailed in Algorithm 1, for two different reasons.
First, the VAA and the $\text{VAA}_{M,\varepsilon}$ are not differentiable
because of the Kronecker delta, which prevents from computing the gradient.
Second, it is likely that hidden states convergence is slow when several RNNs
are stacked. Indeed, the first layers must have reached stability for the
following one to receive a stable input.

In order to solve the first problem, we introduce a differentiable proxy
$\text{VAA}^{*}_{M,\varepsilon}$ of the $\text{VAA}_{M,\varepsilon}$. Instead
of the denominator

$C_{i,j}=\delta\left(\left\lVert
f^{M}(\mathbf{x}_{i},\mathbf{u})-f^{M}(\mathbf{x}_{j},\mathbf{u})\right\lVert\leq\varepsilon\right),$
(6)

that is equal to $1$ when the final states after truncated convergence are
close enough, we use

$C_{i,j}^{*}=1-\frac{\max(0,\left\lVert\tanh
f^{M}(\mathbf{x}_{i},\mathbf{u})-\tanh
f^{M}(\mathbf{x}_{j},\mathbf{u})\right\lVert-\varepsilon)}{\left\lVert\tanh
f^{M}(\mathbf{x}_{i},\mathbf{u})-\tanh
f^{M}(\mathbf{x}_{j},\mathbf{u})\right\lVert}.$ (7)

We note that the value of $C^{*}_{i,j}$ is strictly equal to $1$ if
$f^{M}(\mathbf{x}_{i},u)$ is close enough in Euclidian distance to
$f^{M}(\mathbf{x}_{j},u)$. On the other hand, $C^{*}_{i,j}$ will be close to
$0$ when they are far away. We also note that $C^{*}_{i,j}$ will never be
strictly equal to $0$, but will get closer as the distance increases, since
the fraction tends towards $1$. It can be noted that we are not interested in
states being far apart from each other, but just in them being different.
However, we noticed in the experiments that this small bias provides a good
direction for the gradient in order to reach multistability. For this same
reason, we need to apply a saturating function (hyperbolic tangent in this
case) to the states in order to avoid extreme states when maximizing VAA*. The
resulting differentiable proxy of the VAA is given by

$\text{VAA}^{*}_{M,\varepsilon}(f,\mathcal{X},\mathbf{u})=\frac{1}{\left\lvert\mathcal{X}\right\rvert}\sum^{\left\lvert\mathcal{X}\right\rvert}_{i=1}\dfrac{1}{\overset{\left\lvert\mathcal{X}\right\rvert}{\underset{j=1}{\sum}}\,1-\frac{\max(0,\left\lVert\tanh
f^{M}(\mathbf{x}_{i},\mathbf{u})-\tanh
f^{M}(\mathbf{x}_{j},\mathbf{u})\right\lVert-\varepsilon)}{\left\lVert\tanh
f^{M}(\mathbf{x}_{i},\mathbf{u})-\tanh
f^{M}(\mathbf{x}_{j},\mathbf{u})\right\lVert}}.$ (8)

For maximizing the number of reachable attractors for a given dataset of input
sequences, we thus propose to maximise the proportion of reachable attractors
by SGD on batches of sequences that are sampled from the dataset, using the
VAA* on random hidden states sampled in each of the sequences of the batch.
More precisely, we use SGD to get the VAA* of each layer as close as possible
to $k\in[0,1]$. In practice, we use $k=0.95$ as this proved empirically to
maximise the number of attractors while avoiding too extreme states that could
arise from the approximation of the VAA with $C^{*}$. The loss used is thus
given by

$\mathcal{L}(\mathbf{v},k)=\frac{1}{L}\sum_{l=1}^{L}(v_{l}-k)^{2}\quad$ (9)

where $v_{l}=\text{VAA}^{*}_{M,\varepsilon}(f_{l},\mathcal{X},\mathbf{u})$ is
the estimated multistability of layer $l$ and $L$ is the number of layers in
the RNN. Maximizing the VAA* of each layer separately allows one to tackle the
problem of layer convergence as identified above. To avoid over-fitting
problems, $M$ is sampled uniformly in $\\{1,\dots,M_{\text{max}}(s)\\}$ at
gradient step $s$, where $M_{\text{max}}(s)=\min(M^{*},1+c\cdot s)$ with
$M^{*}$ the maximum stabilization period and $c$ the stabilization period
increment. This progressive increase is required for reaching multistability
smoothly, avoiding gradients problems. For the supervised learning tasks, the
batches of input sequences are sampled in the training set. For the
reinforcement learning tasks, batches of input sequences are sampled from the
exploration policy. Algorithm 3 details the whole warmup procedure for a
dataset $\mathcal{D}$ of input sequences.

1

2

Parameters : $S\in\mathbb{N}$ the number of gradient steps.

$n\in\mathbb{N}$ the batch size.

$\alpha\in\mathbb{R}^{+}$ the learning rate.

$k\in[0,1]$ the target average $\text{VAA}^{*}_{M,\varepsilon}$.

$M^{*}\in\mathbb{N}$ the maximum stabilization period.

$c\in\mathbb{N}$ the stabilization period increment.

$\varepsilon\in\mathbb{R}^{+}$ tolerance when considering state similarity.

$\theta\in\mathbb{R}^{d_{\theta}}$ the parameters of the network.

$L$ the number of layers in the RNN.

3

Data:
$\mathcal{D}=\\{\mathbf{u}_{1:T_{1}}^{1},\dots,\mathbf{u}_{1:T_{N}}^{N}\\}$ a
training set of $N$ input sequences.

4

5for  _$s=1,\dots,S$_ do

6 Sample a batch $\mathcal{B}$ of $n$ sequences in $\mathcal{D}$ without
replacement $\mathcal{B}\sim\mathcal{U}^{n}(\mathcal{D})$.

7

8 Sample a random hidden state in each input sequence
$\mathcal{X}\leftarrow\textnormal{{RandomHiddenStates}}(\mathcal{B},\theta)$.

9

10 Sample $M\sim\mathcal{U}(\\{1,\dots,\min(M,1+s\cdot c)\\})$.

11

12 for  _$l=1,\dots,L$_ do

13 Sample $\mathbf{u}\sim\mathcal{N}(\mathbf{0},\mathbf{1})$.

14 Set $v_{l}=\text{VAA}^{*}_{M,\varepsilon}(f_{l},\mathcal{X},\mathbf{u})$
where $f_{l}$ is the update function of the $l$th RNN layer.

15

16 Compute loss $L\leftarrow\mathcal{L}(\mathbf{v},k)$ where
$\mathbf{v}=\begin{pmatrix}v_{1}&\cdots&v_{L}\end{pmatrix}$.

17

18 Compute gradient $g\leftarrow\nabla_{\theta}L$ with BPTT (over
stabilization period and input sequence).

19

20 Update parameters $\theta\leftarrow\theta-\alpha g$.

21

22 Update maximum stabilization period $M^{*}\leftarrow M^{*}+c$.

23

Algorithm 3 Warming up an RNN

We show in Figure 8 that the warmup procedure effectively increases the VAA*
of each layer in an RNN. Furthermore, we can also see on the right in Figure 8
that as the warmup procedure is carried out, the true VAA measure of the RNN
is increasing as well, even reaching $1$ as the warmup procedure ends.

Figure 8: Evolution of the VAA* for a two-layer RNN (left and middle) and of
the VAA of the network (right) during warmup. The stabilization period for the
VAA is chosen to $M=$2000$$. This network is warmed up on the denoising
dataset and results were averaged over three runs.

###  5.2 Experiments

To demonstrate the impact of warming up RNNs on information restitution tasks,
sequence classification tasks, and in partially observable RL environment, we
tackle all benchmarks introduced in Section 3. We train the LSTM, GRU and MGU
cells with and without warmup and show that their performance is greatly
improved with warmup. As chrono-initialised LSTMs are known to work well, we
also compare our results to such cells, with and without warmup. The
parameters for warming up can vary depending on architectures and needs, but
we found $\alpha=1e^{-2}$, $c=10$, $S=100$, $n=200$ and $M^{*}=200$ to be a
good choice. For the copy first input benchmark, networks are made up of one
$128$ neurons recurrent layer. For the other benchmarks, networks are made up
of two recurrent layers, each of $256$ neurons. For the information
restitution benchmarks, training sets comprise $40\,000$ samples and
performances are evaluated on test sets generated with $50\,000$ samples. For
sequence classification benchmarks, the training set and test set are chosen
to the usual MNIST training set and test set. All averages and standard
deviations reported were computed over three different seeds and training was
done with the Adam optimizer [Kingma and Ba, 2014] and a learning rate of
$1e^{-3}$.

(a) $T=50$

(b) $T=300$

(c) $T=600$

Figure 9: Mean square error ($\pm$ standard deviation) of different
architectures with and without warmup on the test set for the copy first input
benchmark. Results are shown after $50$ epochs and for different values of
$T$.

#### Copy first input benchmark

As can be seen from Figure 9, warming up RNNs greatly improves performances in
the copy first input benchmark, for any sequence length
$T\in\\{50,300,600\\}$. Indeed, classically initialised RNNs have an average
loss above $0.500$ after $50$ epochs, while all warmed-up RNNs have an average
loss below $0.001$ after $50$ epochs. It can also be noted that the chrono-
initialised LSTM are not significantly better when warmed up.

(a) $N=5$

(b) $N=100$

Figure 10: Mean square error ($\pm$ standard deviation) of different
architectures with and without on the denoising benchmark’s test set. Results
are shown with and without constraint on the location of relevant inputs and
after $50$ epochs. Relevant inputs cannot appear in the $N$ last timesteps. In
this experiment, results were obtained with $T=200$.

#### Denoising benchmark

As far as the denoising benchmark is concerned, Figure 10 shows that warmed-up
cells always perform better than classically initialised ones, on sequences of
length $T=200$. However, it can be noted that the average loss is still quite
significant after $50$ epochs for the LSTM and MGU cells, in the case of a
forgetting period of $N=100$. As for the copy first input benchmark, the
chrono-initialised cells seem to perform worse than the other cells, even when
warmed-up. It suggests that the chrono initialisation interacts
disadvantageously with the warmup procedure.

Figure 11: Evolution of the mean cumulative reward obtained by warmed-up and
non-warmed-up agents during their training (left) and mean number of episodes
required to reach the optimal policy (right) on T-Mazes of length 20 (top),
100 (middle) and 200 (bottom).

#### T-Maze benchmark

On the left in Figure 11, we can see the evolution of the expected cumulative
reward of the DRQN policy for the T-Maze environment as a function of the
number of episodes of interaction. It is more than clear that all warmed-up
cells and bistable cells (i.e., BRC and nBRC), are better than the classically
initialised ones on this RL benchmark. As for the other benchmarks, the
chrono-initialised LSTMs seem to interact disadvantageously with the warmup
procedure. In any case, it can be noted that the chrono-initialised LSTMs are
always among the worse cells for this benchmark, with and without warmup.
Furthermore, we can see that warming up cells improves their performance even
more as the length of the T-Maze increases, suggesting that the warmup
procedure and the multistability of an RNN indeed help to tackle tasks with
long-term dependencies. On the right in Figure 11, we can see the number of
episodes required to reach the optimal policy for each cell. It is clear that
warming up a cell speeds up the convergence towards the optimal policy when
time dependencies become large. Indeed, for $L=200$, all warmed-up cells reach
the optimal policy before any classically initialised cell, except for the
chrono-initialised LSTM.

Figure 12: Overall accuracy ($\pm$ standard deviation) on the permuted
sequential MNIST benchmark’s test set after $70$ epochs for different warmup
methods and for different cell types.

#### Permuted sequential MNIST

In Figure 12, we can see the test accuracies after $70$ epochs on the permuted
sequential MNIST benchmark. It is clear that the warmup initialisation does
not help in this task. For the LSTM and GRU, the warmed-up cells are even
worse than the classic cells. This confirms that some tasks such as this
sequence classification benchmark needs more transient dynamics instead of
multistable ones.

(a) $N=72$

(b) $N=472$

Figure 13: Overall accuracy ($\pm$ standard deviation) on permuted sequential-
line MNIST test set after $70$ epochs for different cells with or without
warmup. Images are fed to the recurrent network line by line and $N$ black
lines are added at the bottom of the image after permutation. We note that
when $N$ equals $72$ ($472$) the resulting image has $100$ ($500$) lines.

#### Permuted sequential-line MNIST

In Figure 13, we can see the accuracies of each cell after $70$ epochs on the
test set of the permuted sequential-line MNIST benchmark. For a sequence
length of $100$ (i.e., $N=72$), it is clear that the classically initialised
cells are better at this task. As for the permuted sequential MNIST, this
shows that transient dynamics are important for those sequence classification
tasks, as opposed to information restitution tasks.

###  5.3 Recurrent double-layers

As shown in the previous section and mentioned in the literature [Sussillo and
Barak, 2013], the importance of the transient dynamics of RNNs should not be
neglected for prediction. Indeed, it is easy to see why transient dynamics can
be of importance when trying to tackle a regression task. If information is
only stored in the form of attractors, then there can only be a limited number
of states the network can take, making it very hard to get precise
predictions. We observe that when warming up neural networks they tend to lose
predictive accuracy, at the benefit of easier training on longer sequences.
This leads one to think that RNNs should be built to have both rich transient
and multistability dynamics. We thus propose using a double-layer architecture
that allows one to get precise predictions while maintaining the benefits of
warmup. We simply split each recurrent layer in two equal parts and only
warmup one of them. This allows to endow some part of each layer with
multistability, while the other remains mono-stable with richer transient
dynamics. A double-layer structure is depicted in Figure 14.

InputRNNWarmed-up RNN$\oplus$RNNWarmed-up RNN$\dots$Layer 1Layer 2

Figure 14: Double layer architecture

(a) $T=50$

(b) $T=300$

(c) $T=600$

Figure 15: Mean square error ($\pm$ standard deviation) of different
architectures and different warmup strategies on the test set for the copy
first input benchmark. Results are shown after $50$ epochs and for different
values of $T$.

(a) $N=5$

(b) $N=100$

Figure 16: Mean square error ($\pm$ standard deviation) of different
architectures on the denoising benchmark’s test set. Results are shown with
and without constraint on the location of relevant inputs and after $50$
epochs. Relevant inputs cannot appear in the $N$ last timesteps. In this
experiment, results were obtained with $T=200$. Figure 17: Overall accuracy
($\pm$ standard deviation) on the permuted sequential MNIST benchmark’s test
set after $70$ epochs for different warmup methods and for different cell
types.

(a) $N=72$

(b) $N=472$

Figure 18: Overall accuracy ($\pm$ standard deviation) on permuted sequential-
line MNIST test set after $70$ epochs for different cells with or without
warmup. Images are fed to the recurrent network line by line and $N$ black
lines are added at the bottom of the image after permutation. We note that
when $N$ equals $72$ ($472$) the resulting image has $100$ ($500$) lines.

Figure 19: Evolution of the validation loss on the denoising benchmark for GRU
(left) and LSTM (right). For each cell, three versions have been tested: the
classical one, the warmed-up one and the double-layer one.

As can be seen from Figure 15, Figure 16, Figure 17 and Figure 18, the double-
layer architecture is always among the best performing architecture, for all
four supervised learning benchmarks and for the LSTM, GRU and MGU cells. Even
the chrono-initialised LSTMs are better with the double-layer architecture on
all benchmarks except the copy first input benchmark. It shows that the
double-layer architecture combines both the transient and multistable features
of an RNN. In addition, we can see in Figure 18 that the double-layer
architecture is significantly better than the other architecture, for all
types of cell, on the permuted sequential-line MNIST benchmark with a
forgetting length of $N=472$, a problem that requires both transient and
multistable dynamics.

In Figure 19, we can visualise the evolution of the validation loss on the
denoising benchmark for the GRU and LSTM cells, with the three architectures
(i.e., classic, warmed up and double). It is clear that the warmed-up and
double-layer architectures are better. Additionally, we can see that the
double-layer architecture is significantly faster at learning this task.

##  6 Conclusion

In this work, we introduced a new initialisation procedure, called warmup,
that improve the ability of recurrent neural networks to learn long-term
dependencies. This procedure is motivated by recent work that showed the
importance of fixed points and attractors for the prediction process of
trained RNNs. More precisely, we introduced a lightweight measure called VAA,
that can be optimised at initialisation to endow RNNs with multistable
dynamics. Warmup can be used with any type of recurrent cell and we show that
it vastly improves their performance on problems with long-term dependencies.
In addition, we introduced a new architecture that combines transient and
multistable dynamics through partial warmup. This architecture was shown to
reach a better performance than both classic and warmed-up cells on several
tasks, including information restitution and sequence classifications tasks.

This work also motivates several future works. First, it can be noted that the
double-layer architecture might be worth exploring with different types of
cell. We showed here that there are benefits of using different types of
initialisation for the same type of cell. This might hint at the possibility
of having similar benefits when combining different types of cell that have
different dynamical properties in a single recurrent neural network.
Furthermore, in this paper we have aimed at maximizing the number of
attractors through warmup before training. We noticed however that in some
rare cases, networks loose multistability properties when training. Using VAA
as a regularization loss to avoid this could be interesting. For online
reinforcement learning too, a regularization loss throughout the learning
procedure might make more sense than warming up a priori on random
trajectories. Moreover, we note that not all benchmarks would benefit from
warming up. In fact, it is likely that for several benchmarks, having only a
few attractors could be better. In this regard, it would be interesting to try
to warm up in order to reach a specific number of attractors, rather than for
maximizing them. Finally, the warmup procedure maximises reachable
multistability for a particular dataset of input sequences. Warming up on
totally random input sequences would result in a simpler procedure that might
still provide a good initialisation for reaching multistability.

## Acknowledgements

Gaspard Lambrechts and Nicolas Vecoven gratefully acknowledge the financial
support of the _Wallonia-Brussels Federation_ for their FRIA grant. Florent De
Geeter and Gaspard Lambrechts gratefully acknowledge the financial support of
the _Walloon Region_ for Grant No. 2010235 – ARIAC by DW4AI. Computational
resources have been provided by the _Consortium des Équipements de Calcul
Intensif_ (CÉCI), funded by the _Fonds de la Recherche Scientifique de
Belgique_ (F.R.S.-FNRS) under Grant No. 2502011 and by the Walloon Region.

## References

  * Sussillo and Barak [2013] David Sussillo and Omri Barak.  Opening the black box: low-dimensional dynamics in high-dimensional recurrent neural networks.  _Neural computation_ , 25(3):626–649, 2013.
  * Ceni et al. [2020] Andrea Ceni, Peter Ashwin, and Lorenzo Livi.  Interpreting recurrent neural networks behaviour via excitable network attractors.  _Cognitive Computation_ , 12(2):330–356, 2020\.
  * Maheswaranathan et al. [2019] Niru Maheswaranathan, Alex Williams, Matthew Golub, Surya Ganguli, and David Sussillo.  Reverse engineering recurrent networks for sentiment classification reveals line attractor dynamics.  _Advances in neural information processing systems_ , 32, 2019.
  * Katz and Reggia [2017] Garrett E Katz and James A Reggia.  Using directional fibers to locate fixed points of recurrent neural networks.  _IEEE transactions on neural networks and learning systems_ , 29(8):3636–3646, 2017.
  * Pascanu et al. [2013] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio.  On the difficulty of training recurrent neural networks.  In _International conference on machine learning_ , pages 1310–1318. PMLR, 2013.
  * Cho et al. [2014] Kyunghyun Cho, Bart Van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio.  On the properties of neural machine translation: Encoder-decoder approaches.  _arXiv preprint arXiv:1409.1259_ , 2014.
  * Hochreiter and Schmidhuber [1997] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Jing et al. [2019] Li Jing, Caglar Gulcehre, John Peurifoy, Yichen Shen, Max Tegmark, Marin Soljacic, and Yoshua Bengio.  Gated orthogonal recurrent units: On learning to forget.  _Neural computation_ , 31(4):765–783, 2019.
  * Van Der Westhuizen and Lasenby [2018] Jos Van Der Westhuizen and Joan Lasenby.  The unreasonable effectiveness of the forget gate.  _arXiv preprint arXiv:1804.04849_ , 2018.
  * Marichal et al. [2009] Roberto Marichal, JD PiÑeiro, EJ González, and JM Torres.  New approach of recurrent neural network weight initialization.  In _Advances in Computational Algorithms and Data Analysis_ , pages 537–548. Springer, 2009.
  * Vecoven et al. [2021] Nicolas Vecoven, Damien Ernst, and Guillaume Drion.  A bio-inspired bistable recurrent cell allows for long-lasting memory.  _Plos one_ , 16(6):e0252676, 2021.
  * Werbos [1990] Paul J Werbos.  Backpropagation through time: what it does and how to do it.  _Proceedings of the IEEE_ , 78(10):1550–1560, 1990.
  * Chung et al. [2014] Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio.  Empirical evaluation of gated recurrent neural networks on sequence modeling.  _arXiv preprint arXiv:1412.3555_ , 2014.
  * Zhou et al. [2016] Guo-Bing Zhou, Jianxin Wu, Chen-Lin Zhang, and Zhi-Hua Zhou.  Minimal gated unit for recurrent neural networks.  _International Journal of Automation and Computing_ , 13(3):226–234, 2016.
  * Tallec and Ollivier [2018] Corentin Tallec and Yann Ollivier.  Can recurrent neural networks warp time?  _arXiv preprint arXiv:1804.11188_ , 2018.
  * Bengio et al. [1993] Yoshua Bengio, Paolo Frasconi, and Patrice Simard.  The problem of learning long-term dependencies in recurrent networks.  In _IEEE international conference on neural networks_ , pages 1183–1188. IEEE, 1993.
  * Doya [1993] Kenji Doya.  Bifurcations of recurrent neural networks in gradient descent learning.  _IEEE Transactions on neural networks_ , 1(75):218, 1993.
  * Bakker [2001] Bram Bakker.  Reinforcement learning with long short-term memory.  _Advances in neural information processing systems_ , 14, 2001.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Smallwood and Sondik [1973] Richard D Smallwood and Edward J Sondik.  The optimal control of partially observable markov processes over a finite horizon.  _Operations research_ , 21(5):1071–1088, 1973\.
  * Kaelbling et al. [1998] Leslie Pack Kaelbling, Michael L Littman, and Anthony R Cassandra.  Planning and acting in partially observable stochastic domains.  _Artificial intelligence_ , 101(1-2):99–134, 1998\.
  * Porta et al. [2004] Josep M Porta, Matthijs TJ Spaan, and Nikos Vlassis.  Value iteration for continuous-state POMDPs, 2004.
  * Hausknecht and Stone [2015] Matthew Hausknecht and Peter Stone.  Deep recurrent q-learning for partially observable mdps.  _arXiv preprint arXiv:1507.06527_ , 2015.

##  Appendix A Partially observable Markov decision process

Formally, a POMDP $P$ is an $8$-tuple
$P=(\mathcal{S},\mathcal{A},\mathcal{O},p_{0},T,R,O,\gamma)$ where
$\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space, and
$\mathcal{O}$ is the observation space. The initial state distribution $p_{0}$
gives the probability $p_{0}(\mathbf{s}_{0})$ of
$\mathbf{s}_{0}\in\mathcal{S}$ being the initial state of the decision
process. The dynamics are described by the transition distribution $T$ that
gives the probability $T(\mathbf{s}_{t+1}\mid\mathbf{s}_{t},\mathbf{a}_{t})$
of $\mathbf{s}_{t+1}\in\mathcal{S}$ being the state resulting from action
$\mathbf{a}_{t}\in\mathcal{A}$ in state $\mathbf{s}_{t}\in\mathcal{S}$. The
reward function $R$ gives the immediate reward
$r_{t}=R(\mathbf{s}_{t},\mathbf{a}_{t},\mathbf{s}_{t+1})$ obtained after each
transition. The observation distribution $O$ gives the probability
$O(\mathbf{o}_{t}\mid\mathbf{s}_{t})$ to get observation
$\mathbf{o}_{t}\in\mathcal{O}$ in state $\mathbf{s}_{t}\in\mathcal{S}$.
Finally, the discount factor $\gamma\in[0,1[$ gives the relative importance of
future rewards.

Taking a sequence of $t$ actions ($\mathbf{a}_{0:t-1}$) in the POMDP
conditions its execution and provides a sequence of $t+1$ observations
($\mathbf{o}_{0:t}$). Together, they compose the history
$\eta_{0:t}=(\mathbf{o}_{0:t},\mathbf{a}_{0:t-1})\in\mathcal{H}_{0:t}$ until
timestep $t$, where $\mathcal{H}_{0:t}$ is the set of such histories. Let
$\eta\in\mathcal{H}$ denote a history of arbitrary length sampled in the
POMDP, and let $\mathcal{H}=\bigcup_{t=0}^{\infty}\mathcal{H}_{0:t}$ denote
the set of histories of arbitrary length.

A policy $\pi\in\Pi$ in a POMDP is a mapping from histories to actions, where
$\Pi=\mathcal{H}\rightarrow\mathcal{A}$ is the set of such mappings. A policy
$\pi^{*}\in\Pi$ is said to be optimal when it maximises the expected
discounted sum of future rewards starting from any history
$\eta_{0:t}\in\mathcal{H}_{0:t}$ at time $t\in\mathbb{N}_{0}$

$\pi^{*}\in\operatorname*{arg\,max}_{\pi\in\Pi}\mathop{\underset{\begin{subarray}{c}\pi,P\end{subarray}}{\mathbb{E}}\left[\sum_{t^{\prime}=t}^{\infty}\gamma^{t^{\prime}-t}r_{t^{\prime}}\;\middle|\;\eta_{0:t}\right]},\;\forall\eta_{0:t}\in\mathcal{H}_{0:t},\;\forall
t\in\mathbb{N}_{0}.$ (10)

The history-action value function, or $\mathcal{Q}$-function, is defined as
the maximal expected discounted reward that can be gathered, starting from a
history $\eta_{0:t}\in\mathcal{H}_{0:t}$ at time $t\in\mathbb{N}_{0}$ and an
action $\mathbf{a}_{t}\in\mathcal{A}$

$\mathcal{Q}(\eta_{0:t},\mathbf{a}_{t})=\max_{\pi\in\Pi}\mathop{\underset{\begin{subarray}{c}\pi,P\end{subarray}}{\mathbb{E}}\left[\sum_{t^{\prime}=t}^{\infty}\gamma^{t^{\prime}-t}r_{t^{\prime}}\;\middle|\;\eta_{0:t},\mathbf{a}_{t}\right]},\;\forall\eta_{0:t}\in\mathcal{H}_{0:t},\;\forall\mathbf{a}_{t}\in\mathcal{A},\;\forall
t\in\mathbb{N}_{0}.$ (11)

The $\mathcal{Q}$-function is also the unique solution of the Bellman equation
[Smallwood and Sondik, 1973, Kaelbling et al., 1998, Porta et al., 2004]

$\mathcal{Q}(\eta,\mathbf{a})=\mathop{\underset{\begin{subarray}{c}P\end{subarray}}{\mathbb{E}}\left[r+\gamma\max_{\mathbf{a}^{\prime}\in\mathcal{A}}\mathcal{Q}(\eta^{\prime},\mathbf{a}^{\prime})\;\middle|\;\eta,\mathbf{a}\right]},\;\forall\eta\in\mathcal{H},\;\forall\mathbf{a}\in\mathcal{A}$
(12)

where $\eta^{\prime}=\eta\cup(\mathbf{a},\mathbf{o}^{\prime})$ and $r$ is the
immediate reward obtained when taking action $\mathbf{a}$ in history $\eta$.
From (10) and (11), it can be noticed that any optimal policy satisfies

$\pi^{*}(\eta)\in\operatorname*{arg\,max}_{\mathbf{a}\in\mathcal{A}}\mathcal{Q}(\eta,\mathbf{a}),\;\forall\eta\in\mathcal{H}.$
(13)

##  Appendix B Deep Recurrent Q-learning

The DRQN [Hausknecht and Stone, 2015] algorithm aims at learning a parametric
approximation $\mathcal{Q}_{\theta}$ of the $\mathcal{Q}$-function, where
$\theta\in\mathbb{R}^{d_{\theta}}$ is the parameter vector of a recurrent
neural network. This algorithm is motivated by equation (13) that shows that
an optimal policy can be derived from the $\mathcal{Q}$-function. The strategy
consists of minimising with respect to $\theta$, for all $(\eta,\mathbf{a})$,
the distance between the estimation $\mathcal{Q}_{\theta}(\eta,\mathbf{a})$ of
the LHS of equation (12), and the estimation of the expectation
$\mathbb{E}_{P}[r+\gamma\max_{\mathbf{a}^{\prime}\in\mathcal{A}}\mathcal{Q}_{\theta}(\eta^{\prime},\mathbf{a}^{\prime})]$
of the RHS of equation (12). This is done by using transitions
$(\eta,\mathbf{a},r,\mathbf{o}^{\prime},\eta^{\prime})$ sampled in the POMDP,
with $\eta^{\prime}=\eta\cup(\mathbf{a},\mathbf{o}^{\prime})$.

In practice, this algorithm interleaves the generation of episodes and the
update of the estimation $\mathcal{Q}_{\theta}$. Indeed, in the DRQN
algorithm, the episodes are generated with the $\varepsilon$-greedy policy
derived from the current estimation $\mathcal{Q}_{\theta}$. This stochastic
policy selects actions according to
$\operatorname*{arg\,max}_{\mathbf{a}\in\mathcal{A}}\mathcal{Q}_{\theta}(\cdot,\mathbf{a})$
with probability $1-\varepsilon$, and according to an exploration policy with
probability $\varepsilon$. This exploration policy is defined by a probability
distribution $\mathcal{E}(\mathcal{A})\in\mathcal{P}(\mathcal{A})$ over the
actions, where $\mathcal{P}(\mathcal{A})$ is the set of probability measures
over the action space $\mathcal{A}$. The DRQN algorithm also introduces a
truncation horizon $H$ such that the histories generated in the POMDP have a
maximum length of $H$. Moreover, a replay buffer of histories is used and the
gradient is evaluated on a batch of histories sampled from this buffer.
Furthermore, the parameters $\theta$ are updated with the Adam algorithm
[Kingma and Ba, 2014]. Finally, the target
$r_{t}+\gamma\max_{\mathbf{a}\in\mathcal{A}}\mathcal{Q}_{\theta^{\prime}}(\eta_{0:t+1},\mathbf{a})$
is computed using a past version $\mathcal{Q}_{\theta^{\prime}}$ of the
estimation $\mathcal{Q}_{\theta}$ with parameters $\theta^{\prime}$ that are
updated to $\theta$ less frequently, which eases the convergence towards the
target, and ultimately towards the $\mathcal{Q}$-function.

The DRQN training procedure is detailed in Algorithm 4. In this algorithm, the
output of the RNN is
$\mathbf{y}_{t}=g(\mathbf{h}_{t};\theta)\in\mathbb{R}^{|\mathcal{A}|}$, and it
gives
$\mathcal{Q}_{\theta}(\eta_{0:t},\mathbf{a}),\;\forall\mathbf{a}\in\mathcal{A}$.
The hidden states are given by
$\mathbf{h}_{k}=f(\mathbf{h}_{k-1},\mathbf{x}_{k};\theta),\;\forall
k\in\mathbb{N}_{0}$, with the inputs given by
$\mathbf{x}_{k}=(\mathbf{a}_{k-1},\mathbf{o}_{k}),\;\forall t\in\mathbb{N}$
and $\mathbf{x}_{0}=(\mathbf{0},\mathbf{o}_{0})$. From the approximation
$\mathcal{Q}_{\theta}$, the policy $\pi_{\theta}$ is given by
$\pi_{\theta}(\eta)=\operatorname*{arg\,max}_{\mathbf{a}\in\mathcal{A}}\mathcal{Q}_{\theta}(\eta,\mathbf{a})$.

1

2

Parameters : $N\in\mathbb{N}$ the buffer capacity.

$C\in\mathbb{N}$ the target update period in term of episodes.

$E\in\mathbb{N}$ the number of episodes.

$H\in\mathbb{N}$ the truncation horizon.

$I\in\mathbb{N}$ the number of gradient steps after each episode.

$\varepsilon\in\mathbb{R}$ the exploration rate.

$\mathcal{E}(\mathcal{A})\in\mathcal{P}(\mathcal{A})$ the exploration policy
probability distribution.

$\alpha\in\mathbb{R}$ the learning rate.

$B\in\mathbb{N}$ the batch size.

$\theta\in\mathbb{R}^{d_{\theta}}$ the initial parameters of the network.

$\theta^{\prime}\in\mathbb{R}^{d_{\theta}}$ the initial parameters of the
target network.

3

4Initialise weights $\theta$ randomly

5

6Fill replay buffer $\mathcal{B}$ with random transitions from the exploration
policy $\mathcal{E}(\mathcal{A})$.

7 if  _warmup_ then

8 Let $\mathcal{D}$ be the set of histories $\eta$ (input sequences) in replay
buffer $\mathcal{B}$.

9 Warmup ($\mathcal{D},\theta$) using default parameters of the Warmup
algorithm.

10for  _$e=0,\dots,E-1$_ do

11

12 if  _$e\bmod C=0$_ then

13 Update target network with $\theta^{\prime}\leftarrow\theta$

// Generate new episode, store history and rewards

14

15 Draw an initial state $\mathbf{s}_{0}$ according to $p_{0}$ and observe
$\mathbf{o}_{0}$

16

17 Let $\eta_{0:0}=(\mathbf{o}_{0})$

18

19 for  _$t=0,\dots,H-1$_ do

20 Select $\mathbf{a}_{t}\sim\mathcal{E}(\mathcal{A})$ with probability
$\varepsilon$, otherwise select
$\mathbf{a}_{t}=\operatorname*{arg\,max}_{\mathbf{a}\in\mathcal{A}}\left\\{\mathcal{Q}_{\theta}(\eta_{0:t},\mathbf{a})\right\\}$

21

22 Take action $\mathbf{a}_{t}$ and observe $r_{t}$ and $\mathbf{o}_{t+1}$

23 Let
$\eta_{0:t+1}=(\mathbf{o}_{0},\mathbf{a}_{0},\mathbf{o}_{1},\dots,\mathbf{o}_{t+1})$

24

25 if  _$|\mathcal{B}| <N$_ then add
$(\eta_{0:t},\mathbf{a}_{t},r_{t},\mathbf{o}_{t+1},\eta_{0:t+1})$ in replay
buffer $\mathcal{B}$

26 else replace oldest transition in replay buffer $\mathcal{B}$ by
$(\eta_{0:t},\mathbf{a}_{t},r_{t},\mathbf{o}_{t+1},\eta_{0:t+1})$

27

28 if  _$\mathbf{o}_{t+1}$ is terminal_ then

29 break

// Optimise recurrent Q-network

30

31 for  _$i=0,\dots,I-1$_ do

32

33 Sample $B$ transitions
$(\eta_{0:t}^{b},\mathbf{a}_{t}^{b},r_{t}^{b},\mathbf{o}_{t+1}^{b},\eta_{0:t+1}^{b})$
uniformly from the replay buffer $\mathcal{B}$

34

35 Compute targets
$y^{b}=\begin{cases}r_{t}^{b}+\gamma\max_{\mathbf{a}\in\mathcal{A}}\left\\{\mathcal{Q}_{\theta^{\prime}}(\eta_{0:t+1}^{b},\mathbf{a})\right\\}&\text{if
}\mathbf{o}_{t+1}^{b}\text{ is not terminal}\\\
r_{t}^{b}&\text{otherwise}\end{cases}$

36

37 Compute loss
$L=\sum_{b=0}^{B-1}\left(y^{b}-\mathcal{Q}_{\theta}(\eta_{0:t}^{b},\mathbf{a}_{t}^{b})\right)^{2}$

38

39 Compute direction $g$ using Adam optimiser, perform gradient step
$\theta\leftarrow\theta+\alpha g$

40

Algorithm 4 DRQN - $\mathcal{Q}$-function approximation

##  Appendix C T-Maze environment

The T-Maze environment is a POMDP
$(\mathcal{S},\mathcal{A},\mathcal{O},p_{0},T,R,O,\gamma)$ parameterised by
the maze length $L\in\mathbb{N}$. The formal definition of this environment is
given below.

$\mathbf{m}=\text{Up}$$\mathbf{m}=\text{Down}$$\mathbf{c}=(0,0)$$\mathbf{c}=(1,0)$$\mathbf{c}=(2,0)$$\mathbf{c}=(3,0)$$\mathbf{c}=(4,0)$$\mathbf{c}=(5,0)$$\mathbf{c}=(6,0)$$\dots$$\mathbf{c}=(L,0)$$\mathbf{c}=(L,1)$$\mathbf{c}=(L,-1)$$\mathbf{c}=(0,0)$$\mathbf{c}=(1,0)$$\mathbf{c}=(2,0)$$\mathbf{c}=(3,0)$$\mathbf{c}=(4,0)$$\mathbf{c}=(5,0)$$\mathbf{c}=(6,0)$$\dots$$\mathbf{c}=(L,0)$$\mathbf{c}=(L,1)$$\mathbf{c}=(L,-1)$

Figure 20: T-Maze state space. Initial states in blue, terminal states in
grey, and treasure states hatched.

#### State space

The discrete state space $\mathcal{S}$ is composed of the set of positions
$\mathcal{C}$ for the agent in each of the two maze layouts $\mathcal{M}$. The
maze layout determines the position of the treasure. Formally, we have

$\displaystyle\mathcal{S}$ $\displaystyle=\mathcal{M}\times\mathcal{C}$ (14)
$\displaystyle\mathcal{M}$ $\displaystyle=\\{\text{Up},\text{Down}\\}$ (15)
$\displaystyle\mathcal{C}$
$\displaystyle=\left\\{(0,0),\dots,(L,0)\right\\}\cup\\{(L,1),(L,-1)\\}$ (16)

A state $\mathbf{s}_{t}\in\mathcal{S}$ is thus defined by
$\mathbf{s}_{t}=(\mathbf{m}_{t},\mathbf{c}_{t})$ with
$\mathbf{m}_{t}\in\mathcal{M}$ and $\mathbf{c}_{t}\in\mathcal{C}$. Let us also
define
$\mathcal{F}=\left\\{\mathbf{s}_{t}=(\mathbf{m}_{t},\mathbf{c}_{t})\in\mathcal{S}\mid\mathbf{c}_{t}\in\\{(L,1),(L,-1)\\}\right\\}$
the set of terminal states, four in number.

#### Action space

The discrete action space $\mathcal{A}$ is composed of the four possible moves
that the agent can take

$\mathcal{A}=\left\\{(1,0),(0,1),(-1,0),(0,-1)\right\\}$ (17)

that correspond to Right, Up, Left and Down, respectively.

#### Observation space

The discrete observation space $\mathcal{O}$ is composed of the four partial
observations of the state that the agent can perceive

$\mathcal{O}=\left\\{\text{Up},\text{Down},\text{Corridor},\text{Junction}\right\\}.$
(18)

#### Initial state distribution

The two possible initial states are
$\mathbf{s}_{0}^{\text{Up}}=(\text{Up},(0,0))$ and
$\mathbf{s}_{0}^{\text{Down}}=(\text{Down},(0,0))$, depending on the maze in
which the agent lies. The initial state distribution
$p_{0}:\mathcal{S}\rightarrow[0,1]$ is thus given by

$p_{0}(\mathbf{s}_{0})=\begin{cases}0.5&\text{ if
}\mathbf{s}_{0}=\mathbf{s}_{0}^{\text{Up}}\\\ 0.5&\text{ if
}\mathbf{s}_{0}=\mathbf{s}_{0}^{\text{Down}}\\\ 0&\text{ otherwise
}\end{cases}$ (19)

#### Transition distribution

The transition distribution function
$T\colon\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow[0,1]$ is
given by

$T(\mathbf{s}_{t+1}\mid\mathbf{s}_{t},\mathbf{a}_{t})=\delta_{f(\mathbf{s}_{t},\mathbf{a}_{t})}(\mathbf{s}_{t+1})$
(20)

where $\mathbf{s}_{t}\in\mathcal{S},\mathbf{a}_{t}\in\mathcal{A}$ and
$\mathbf{s}_{t+1}\in\mathcal{S}$, and $f$ is given by

$f(\mathbf{s}_{t},\mathbf{a}_{t})=\begin{cases}\mathbf{s}_{t+1}=(\mathbf{m}_{t},\mathbf{c}_{t}+\mathbf{a}_{t})&\text{
if
}\mathbf{s}_{t}\not\in\mathcal{F},\mathbf{c}_{t}+\mathbf{a}_{t}\in\mathcal{C}\\\
\mathbf{s}_{t+1}=(\mathbf{m}_{t},\mathbf{c}_{t})&\text{ otherwise
}\end{cases}$ (21)

where $\mathbf{s}_{t}=(\mathbf{m}_{t},\mathbf{c}_{t})\in\mathcal{S}$ and
$\mathbf{a}_{t}\in\mathcal{A}$.

#### Reward function

The reward function
$R:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}$ is
given by

$R(\mathbf{s}_{t},\mathbf{a}_{t},\mathbf{s}_{t+1})=\begin{cases}0&\text{ if
}\mathbf{s}_{t}\in\mathcal{F}\\\ 0&\text{ if
}\mathbf{s}_{t}\not\in\mathcal{F},\mathbf{s}_{t+1}\not\in\mathcal{F},\mathbf{s}_{t}\neq\mathbf{s}_{t+1}\\\
-0.1&\text{ if
}\mathbf{s}_{t}\not\in\mathcal{F},\mathbf{s}_{t+1}\not\in\mathcal{F},\mathbf{s}_{t}=\mathbf{s}_{t+1}\\\
4&\text{ if
}\mathbf{s}_{t}\not\in\mathcal{F},\mathbf{s}_{t+1}\in\mathcal{F},\mathbf{c}_{t+1}=\begin{cases}(L,1)&\text{
if }\mathbf{m}_{t+1}=\text{Up}\\\ (L,-1)&\text{ if
}\mathbf{m}_{t+1}=\text{Down}\end{cases}\\\ -0.1&\text{ if
}\mathbf{s}_{t}\not\in\mathcal{F},\mathbf{s}_{t+1}\in\mathcal{F},\mathbf{c}_{t+1}=\begin{cases}(L,-1)&\text{
if }\mathbf{m}_{t+1}=\text{Up}\\\ (L,+1)&\text{ if
}\mathbf{m}_{t+1}=\text{Down}\end{cases}\end{cases}$ (22)

where
$\mathbf{s}_{t}=(\mathbf{m}_{t},\mathbf{c}_{t})\in\mathcal{S},\mathbf{a}_{t}\in\mathcal{A}$
and $\mathbf{s}_{t+1}=(\mathbf{m}_{t+1},\mathbf{c}_{t+1})\in\mathcal{S}$.

#### Observation distribution

In the T-Maze, the observations are deterministic. The observation
distribution $O\colon\mathcal{S}\times\mathcal{O}\rightarrow[0,1]$ is given by

$O(\mathbf{o}_{t}\mid\mathbf{s}_{t})=\begin{cases}1&\text{ if
}\mathbf{o}_{t}=\text{Up},\mathbf{c}_{t}=(0,0),\mathbf{m}_{t}=\text{Up}\\\
1&\text{ if
}\mathbf{o}_{t}=\text{Down},\mathbf{c}_{t}=(0,0),\mathbf{m}_{t}=\text{Down}\\\
1&\text{ if
}\mathbf{o}_{t}=\text{Corridor},\mathbf{c}_{t}\in\left\\{(1,0),\dots,(L-1,0)\right\\}\\\
1&\text{ if
}\mathbf{o}_{t}=\text{Junction},\mathbf{c}_{t}\in\left\\{(L,0),(L,1),(L,-1)\right\\}\\\
0&\text{ otherwise }\end{cases}$ (23)

where $\mathbf{s}_{t}=(\mathbf{m}_{t},\mathbf{c}_{t})\in\mathcal{S}$ and
$\mathbf{o}_{t}\in\mathcal{O}$.

#### Exploration policy

The exploration policy $\mathcal{E}:\mathcal{A}\rightarrow[0,1]$ is a
stochastic policy that is given by $\mathcal{E}(\text{Right})=1/2$ and
$\mathcal{E}(\text{Other})=1/6$ where
$\text{Other}\in\left\\{\text{Up},\text{Left},\text{Down}\right\\}$. It
enforces to explore the right of the maze layouts. This exploration policy,
tailored to the T-Maze environment, allows one to speed up the training
procedure, without interfering with the study of this work.

#### Truncation horizon

The truncation horizon $H$ of the DRQN algorithm is chosen such that the
expected displacement of an agent moving according to the exploration policy
in a T-Maze with an infinite corridor on both sides is greater than $L$. Let
$r=\mathcal{E}(\text{Right})$ and $l=\mathcal{E}(\text{Left})$. In this
infinite T-Maze, starting at $0$, the expected displacement after one timestep
is $\bar{x}_{1}=r-l$. By independence, $\bar{x}_{H}=H\bar{x}_{1}$ such that,
for $\bar{x}_{H}\geq L$, the time horizon is given by

$H=\left\lceil\frac{L}{r-l}\right\rceil.$ (24)


# The Higgs boson as a self-similar system: A new solution to the hierarchy
problem

M. Ahmadvand<EMAIL_ADDRESS>School of Particles and Accelerators, Institute
Research in Fundamental Sciences (IPM), P. O. Box 19395-5531, Tehran, Iran

###### Abstract

We propose a new solution to the hierarchy (naturalness) problem, concerning
quantum corrections of the Higgs mass. Suggesting the Higgs boson as a system
with a self-similar internal structure, we calculate its two-point function
and find that the quadratic divergence is replaced by a logarithmic one. It is
shown that the partonic-like distribution follows the Tsallis statistics and
also high energy physics experimental data for the Higgs transverse momentum
distribution can be described by the self-similar statistical model.

Introduction.— After the discovery of the Higgs boson with a mass around
$125\,\mathrm{GeV}$ at the Large Hadron Collider (LHC), although the last
piece of the Standard Model (SM) has been found ATLAS:2012yve ; CMS:2012qbp ;
CMS:2013btf , some features of the Higgs boson are still under debate and
investigation. Moreover, there are important issues, such as the mass of
neutrinos, matter-antimatter asymmetry, dark matter and the hierarchy problem,
which are left unanswered in the SM and should be addressed in a more
fundamental theory. On the other hand, no significant deviations from the SM
predictions have been observed so far at high energy collisions and the theory
may be generalized to high energy scales. In this sense, the hierarchy problem
of scales between the weak and higher energy scales is more evident when
considering quantum fluctuations to calculate the Higgs squared mass
corrections Veltman:1980mj

$\delta
m_{h}^{2}\sim\frac{\Lambda^{2}}{v^{2}}\left(m_{h}^{2}+2m_{W}^{2}+m_{Z}^{2}-4m_{t}^{2}\right),$
(1)

where $\Lambda$ is the cutoff, the highest accessible energy scale, $v\simeq
246\,\mathrm{GeV}$ is the electroweak (EW) symmetry breaking scale and
$m_{h}$, $m_{W,Z}$, and $m_{t}$ denote the mass of Higgs, gauge boson and top
quark particles. Thus, if $\Lambda$ is very large, for instance as larger as
the Planck mass, the corrections will be extremely greater than the Higgs mass
value.

Theories beyond the SM may have an additional contribution to the Higgs mass
such that it is finely adjusted to cancel $\delta m_{h}^{2}$. Supersymmetric
models Draper:2016pys  and composite Higgs scenarios Panico:2015jxa  are such
endeavors to avoid this problem. The fine-tuning in cancellation can be
measured as

$\Delta\geq\frac{\delta
m_{h}^{2}}{m_{h}^{2}}\sim\left(\frac{\Lambda}{450\,\mathrm{GeV}}\right)^{2}.$
(2)

Therefore, in these new theories, for $\Lambda\gg\mathrm{TeV}$, the Higgs
boson mass cannot be found due to the cancellation.

In the SM, the Higgs boson is the only elementary scalar field and all other
scalars are bound states of the strongly coupled sector. If the Higgs is also
elaborated to be a composite bound state, it can be originated from a new
strongly interacting dynamics, since Quantum Chromodynamics (QCD) cannot be
responsible for its construction. That is the notion established in composite
Higgs models Ahmadvand:2020izy , avoiding the hierarchy problem.

In this paper, the Higgs boson is considered as a complex system with an
internal structure which reveals a statistical self-similarity behavior in
that constituents are similar to the main system at a different level of
scale.

Self-similar objects and patterns, known as fractals, have been widely studied
in various areas of science including mathematics, biology and physics, as
widely found in nature (not necessarily as exact fractals) Mandelbrot . In
particle physics, especially in strong interactions, self-similar systems have
been used to model hadrons which are composed of hadrons Hagedorn:1965st ;
Frautschi:1971ij ; Chew:1971xwl . These models were able to describe many
features of the hadronic system consistent with experimental data.
Additionally, in high energy collision experiments of strong interactions,
showing nearly scale invariant properties, the transverse momentum data would
be well described in this context Bediaga:1999hv ; Cleymans:2011in ;
Deppman:2012us ; Marques:2012px  if the Boltzmann-Gibbs statistics were
replaced by the Tsallis statistics Tsallis:1987eu ; Curado:1991jc .

In the present work, we first obtain the two-point function of the Higgs with
the self-similar internal structure, whose partonic-like distribution is
described by the Tsallis statistics. We then show the modification of the
field strength renormalization gives rise to the logarithmic divergence of the
Higgs mass corrections. Note that the Higgs can be resulted from a new QCD-
like dynamics at high energies and hence not only is the hierarchy problem
addressed in this setup, but also the prediction of the Higgs mass is
naturally feasible without restricting the given confining scale to be around
the TeV.

Correlation function.— We now investigate the two-point correlation function
of the Higgs, in accordance with the aforementioned considerations, for the
interacting theory, $\langle\Omega|\phi(x)\phi(y)|\Omega\rangle$. Using a
complete set of intermediate states at the scale level $n$

$\mathbf{1}=\sum_{\psi_{n}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{2E_{p}}|\psi_{n}\rangle\langle\psi_{n}|,$
(3)

the two-point function for $x^{0}>y^{0}$ can be written as

$\langle\Omega|\phi(x)\phi(y)|\Omega\rangle=\sum_{\psi_{n}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{2E_{p}}\langle\Omega|\phi(x)|\psi_{n}\rangle\langle\psi_{n}|\phi(y)|\Omega\rangle,$
(4)

where $E_{p}=\sqrt{p^{2}+m^{2}}$ and $|\psi_{n}\rangle$ is the self-similar
partonic state which can be an eigenstate of 4-momentum $P^{\mu}$. Because of
translational invariance,

$\displaystyle\langle\Omega|\phi(x)|\psi_{n}\rangle$
$\displaystyle=\langle\Omega|\phi(0)|\psi_{n}\rangle\,e^{-ip.x}$
$\displaystyle\langle\psi_{n}|\phi(y)|\Omega\rangle$
$\displaystyle=\langle\psi_{n}|\phi(0)|\Omega\rangle\,e^{ip.y},$ (5)

thus Eq. (4) is expressed as

$\langle\Omega|\phi(x)\phi(y)|\Omega\rangle=\sum_{\psi_{n}}\int\frac{d^{4}p}{(2\pi)^{4}}\frac{e^{-ip.(x-y)}}{p^{2}-m^{2}+i\epsilon}|\langle\Omega|\phi(0)|\psi_{n}\rangle|^{2}.$
(6)

For $y^{0}>x^{0}$ the same procedure can be written to express the time-
ordered product of the function. Furthermore, we can represent the scalar
function as the K$\ddot{\mathrm{a}}$llen-Lehmann representation in terms of a
spectral function

$\langle\Omega|T\phi(x)\phi(y)|\Omega\rangle=\int_{0}^{\infty}\frac{dM^{2}}{(2\pi)}\rho(M^{2})\,D_{F}(x-y;M^{2}),$
(7)

where $D_{F}$ is the Feynman propagator. Rewriting the effective state in
terms of self-similar partonic constituent states

$|\psi_{n}\rangle=\sum_{i}\langle\psi_{i}|\psi_{n}\rangle|\psi_{i}\rangle,$
(8)

we can explore the effect of internal structure. In this case the spectral
function of the main system will be

$\rho(M^{2})\sim\sum_{i}(2\pi)\delta(M^{2}-m_{h}^{2})|\langle\Omega|\phi(0)|\psi_{i}\rangle|^{2}|\langle\psi_{i}|\psi_{n}\rangle|^{2}.$
(9)

As a result, the two-point function becomes

$\langle\Omega|T\phi(x)\phi(y)|\Omega\rangle\sim\sum_{i}\int_{0}^{\infty}dM^{2}\,\delta(M^{2}-m_{h}^{2})\,Z\,|\langle\psi_{i}|\psi_{n}\rangle|^{2}D_{F},$
(10)

where $Z$ is effectively defined as the field strength renormalization and

$\displaystyle|\langle\psi_{i}|\psi_{n}\rangle|^{2}$
$\displaystyle=\sum_{i}\int d\Pi_{i}~{}|\mathcal{M}|^{2}$
$\displaystyle=\sum_{i}\bigg{(}\prod_{i}\int\frac{d^{4}p_{i}}{(2\pi)^{4}}\bigg{)}|\mathcal{M}_{0}|^{2}\,f\,(2\pi)^{4}\delta^{4}\left(\sum
p_{i}-k\right).$ (11)

Here, $\mathcal{M}$ is the probability amplitude of finding the subsystem
$i$th of the system with mass $m_{i}$ as a fraction of $M$, $\mathcal{M}_{0}$
can be thought of as a scattering amplitude of the partonic system to $N$
constituent fields with a different energy level with respect to the system,
and $f$ includes the contribution of substructures and determines the parton
distribution function of the scale invariant parameter in the model. Hence, it
is expected that $|\mathcal{M}_{0}|^{2}$ is proportional to the dimensional
parameter of the system $M^{2}$ which is normalized to the energy scale of the
larger system $E_{k}^{2}$, that is to say, the accessible energy of the
intermediate state which can be probed,111In the perturbative method with
Feynman diagrams, $k$ is the momentum which is integrated over at the loop
level and here as can be seen in the later analysis, the term associated with
$E_{k}$ is included in the modified field strength renormalization. thereby we
approximate it as $|\mathcal{M}_{0}|^{2}\sim M^{2}/E_{k}^{2}$. Now, one should
find how this scale invariant parameter of the system, $x=M/E_{k}=m_{i}/M$, is
distributed to partons through the distribution function.

Based on the Maxwell distribution for the ideal gas with $N$ particles,

$f(K)=\mu^{-\frac{3N}{2}}K^{\frac{3N}{2}-1}e^{-\frac{K}{\mu}},$ (12)

where $K=\mathbf{k}^{2}/(2M)$ and $\mu$ is the normalization factor, we try to
obtain the distribution for a system having an internal energy and hence with
a total energy $K+I=K(1+I/K)$ where $K$ stands for the kinetic energy of the
system and $I$ for the internal energy. For the self-similar system, $I/K$
remains constant as the scale invariant parameter for all subsystems such that
$I/K=I^{(n)}/K^{(n)}=x$. Due to the self-similarity feature, the energy
distribution of internal subsystems should be equal at different levels and
follow the main structure. Thus, taking the contribution of all subsystems
into account, we consider the total energy $U$ of the system as follows

$U=K\bigg{(}1+\frac{I}{K}\bigg{)}^{\alpha}=K\left(1+x\right)^{\alpha},$ (13)

where $\alpha$ denotes the fractal index and will be determined in the
following by the relevant statistics of the system. From Eq. (13), $U\sim
K+I\alpha(1+x)^{\alpha^{\prime}}$, and for this system we can take
$f(U)=f(K)f(I)f(x)$ as

$\displaystyle f(U)$
$\displaystyle=\mu^{-\frac{3N}{2}}U^{\frac{3N}{2}-1}e^{-\frac{U}{\mu}}f(x)$
$\displaystyle=\mu^{-\frac{3N}{2}}K^{\frac{3N}{2}-1}\left(1+x\right)^{-\frac{3N}{2}\alpha-\tilde{\alpha}}e^{-\frac{K}{\mu}\left(1+x\right)^{\alpha}}$
(14)

where in the last step, we used the self-similarity of the internal structure
which can be also indicated from the following integration. Therefore, the
integration over 3-momentum components of Eq. (The Higgs boson as a self-
similar system: A new solution to the hierarchy problem) for the distribution
is written as

$\displaystyle\sum_{i}\bigg{(}\prod_{i}\int\frac{d^{3}p_{i}}{(2\pi)^{3}}e^{\frac{p_{i}^{2}}{2M}}f(I)f(x)\bigg{)}(2\pi)^{3}\delta^{3}\left(\sum\vec{p_{i}}-\mathbf{k}\right)$
$\displaystyle=\int
dK\,\mu^{-\frac{3N}{2}}K^{\frac{3N}{2}-1}\left(1+x\right)^{-\frac{3N}{2}\alpha-\tilde{\alpha}}e^{-\frac{K}{\mu}\left(1+x\right)^{\alpha}}=\left(1+x\right)^{-\tilde{\alpha}}.$
(15)

Notice that by substituting $\tilde{\alpha}=q/(q-1)$ and
$E_{k}=\tilde{\epsilon}/(q-1)$, it can be seen the distribution,
$[1+(q-1)M/\tilde{\epsilon}]^{-q/(q-1)}$, follows the non-extensive Tsallis
statistics, where $q$ is the entropic factor of the system and for the limit
$q\rightarrow 1$ the Boltzmann-Gibbs is reproduced.

As a result, the probability term, Eq. (The Higgs boson as a self-similar
system: A new solution to the hierarchy problem), will be

$\sum_{i}\bigg{(}\prod_{i}\int\frac{d^{4}p_{i}}{(2\pi)^{4}}\bigg{)}|\mathcal{M}_{0}|^{2}f(2\pi)^{4}\delta^{4}\left(\sum
p_{i}-k\right)=x^{2}\left(1+x\right)^{-\tilde{\alpha}}.$ (16)

Eventually, we obtain the Fourier transformation of the two-point function as
follows

$\displaystyle\int
d^{4}x~{}e^{ip(x-y)}\langle\Omega|T\phi(x)\phi(y)|\Omega\rangle$
$\displaystyle\sim\frac{iZ\,x^{2}(1+x)^{\tilde{\alpha}}}{p^{2}-m_{h}^{2}+i\epsilon}\equiv\frac{i\widetilde{Z}}{p^{2}-m_{h}^{2}+i\epsilon}.$
(17)

Consequently, the effect of the internal structure is included in
$\widetilde{Z}\equiv Z\tilde{f}$ where
$\tilde{f}=(1+m_{h}/E_{k})^{\tilde{\alpha}}\,m_{h}^{2}/E_{k}^{2}$.

Analogous to the sum of one-particle irreducible (1PI) diagrams, denoted by
$\mathbb{M}$, for this system the Higgs two-point function can be expressed as
a geometric series of modified diagrams, $\mathbb{\widetilde{M}}$, so that

$\frac{i\widetilde{Z}}{p^{2}-m_{h}^{2}+i\epsilon}\sim\frac{i}{p^{2}-m_{0h}^{2}-\mathbb{\widetilde{M}}}$
(18)

where $\mathbb{\widetilde{M}}=\tilde{f}\mathbb{M}$ and can be attained from
Eq. (18) by expanding the denominator of the right hand side close to the pole

$\left(p^{2}-m_{h}^{2}\right)\bigg{(}1-\frac{d\mathbb{\widetilde{M}}}{dp^{2}}\bigg{|}_{p^{2}=m_{h}^{2}}\bigg{)}+\mathcal{O}\left((p^{2}-m_{h}^{2})^{2}\right),$
(19)

and comparing it to the left hand side of the equation. Thereby, from
$\widetilde{Z}^{-1}\sim 1-\delta\widetilde{Z}$, we can find
$\delta\widetilde{Z}=\tilde{f}\delta
Z=d\mathbb{\widetilde{M}}/dp^{2}=\tilde{f}d\mathbb{M}/dp^{2}$. This result can
be also consistently applied to the Lehmann-Symanzik-Zimmermann (LSZ) formula
Lehmann:1954rq  in that the sum of 1PI insertions in the propagator is equal
to that of amputed scattering diagrams, implying $\mathcal{M}(p\rightarrow
p)=\widetilde{Z}\mathbb{M}=Z\widetilde{\mathbb{M}}$.

To clarify this effect, we calculate $\widetilde{\mathbb{M}}$ at the one-loop
order, for the Feynman diagram with the self coupling interaction; the same
procedure holds for other Higgs interactions. Since $m_{h}^{2}/E^{2}_{k}<1$,
for the sake of simplicity in the calculation, we consider $\tilde{f}\sim
m_{h}^{2}/E^{2}_{k}$. Thus, for the mentioned diagram

$\displaystyle-i\mathbb{\widetilde{M}}$
$\displaystyle=\frac{-i\lambda}{2}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{m_{h}^{2}}{(k^{2}+m_{h}^{2})}\frac{i}{(k^{2}-m_{h}^{2})}$
$\displaystyle=\frac{-i\lambda
m_{h}^{2}}{32\pi^{2}}\int_{0}^{\Lambda^{2}}\frac{du_{E}\,u_{E}}{u_{E}^{2}+m_{h}^{4}}=\frac{-i\lambda
m_{h}^{2}}{32\pi^{2}}\ln\frac{\Lambda^{2}+m_{h}^{2}}{m_{h}^{2}},$ (20)

where $u=k^{2}$ and $u_{E}$ is used for the Euclidean variable. As can be seen
from Eq. (The Higgs boson as a self-similar system: A new solution to the
hierarchy problem), and generally for other Higgs two-point function diagrams,
we deal with the logarithmic divergence which is canceled by $\delta
m_{h}^{2}$.

Discussion.— As already mentioned, the obtained distribution obeys the non-
extensive Tsallis statistics. The non-extensivity can be realized from the
non-additive entropy Tsallis:1987eu  and the factor $q$ is a measure of non-
additivity.

In high energy experiments, such a distribution has been used for describing
the transverse momentum $p_{\mathrm{T}}$ distribution of hadronic systems
whose experimental data in a good agreement can be modeled by the following
fitting function STAR:2006nmo ; ATLAS:2010jvh ; PHENIX:2011rvu ; ALICE:2011gmo
; CMS:2011jlm

$\frac{dN}{Ndp_{\mathrm{T}}}=\tilde{c}p_{\mathrm{T}}\frac{(n-1)(n-2)}{nC\left(nC+m(n-2)\right)}\bigg{(}1+\frac{\sqrt{p^{2}_{\mathrm{T}}+m^{2}}-m}{nC}\bigg{)}^{-n},$
(21)

where $\tilde{c}$ is the normalization, $n$, $C$ are fitting parameters and
$m$ stands for the hadron mass. The number of events for a given cross section
is $N=\sigma\int\mathcal{L}dt$ and $\int\mathcal{L}dt$ is the integrated
luminosity. The fitting parameters can also be identified in terms of Tsallis
parameters $T_{0}$ and $q$ as $n\equiv q/(q-1)$ and $nC\equiv T_{0}/(q-1)$
Cleymans:2011in .

In addition, in high energy collisions, the Higgs transverse momentum
$p_{\mathrm{T}}^{H}$ is a key observable due to which one can study its
properties and the dynamics of the produced system as well as distortions of
its SM predictions. We try to describe the $p_{\mathrm{T}}$ distribution of
the Higgs system by the mentioned statistical model, the Tsallis function,
with its parameters. We consider the combined differential cross section of
$H\rightarrow\gamma\gamma$ (diphoton), $H\rightarrow ZZ^{*}$ and $H\rightarrow
b\bar{b}$ (a bottom quark-antiquark pair) decay channels for
$p_{\mathrm{T}}^{H}$ reported by the CMS collaboration at
$\sqrt{s}=13\,\mathrm{TeV}$ CMS:2018gwt . (Similar measurements reported by
the ATLAS collaboration can be found in ATLAS:2018pgp .) We fit the spectra
for $d\sigma/dp_{\mathrm{T}}^{H}$ to the Tsallis fit function, fixing
$\tilde{c}$ to the measured total cross section $61.1\,\mathrm{pb}$
CMS:2018gwt . As shown in Fig. (1), the experimental data can be well fitted
to the function for $q=1.79\pm 0.13$ and $T_{0}=3.7\pm 0.8\,\mathrm{GeV}$. The
result can be an evidence for the self-similarity feature of the Higgs and
also using fitting parameters, more detailed investigation may constrain the
mass spectrum of this type of particles in future high energy physics
explorations.

Another aspect of Higgs properties can be studied through the running of its
interaction coupling constants by means of the beta function and the Callan-
Symanik equation. In the present setup, taking into account the modified 1PI
diagrams, corrections can be obtained. For instance at the one-loop order for
the self coupling interaction diagram, the four-point function would be finite
and the relevant contribution to the beta function, with regard to divergent
logarithms up to $\lambda^{2}$ order, is associated with two-loop two-point
function diagrams. However, the result is still compatible to the evolution of
$\lambda$ and related bounds Reina:2012fs .

Figure 1: The combined measurement of differential cross section for
$H\rightarrow\gamma\gamma$, $H\rightarrow ZZ^{*}$ and $H\rightarrow b\bar{b}$
decay channels as a function of $p_{\mathrm{T}}^{H}$ is shown as black points
with error bars CMS:2018gwt . Using the MINUIT package James:1994vla , the
spectra are fitted to the Tsallis function for $q=1.79\pm 0.13$ and
$T_{0}=3.7\pm 0.8\,\mathrm{GeV}$.

Acknowledgments.— The author is so grateful to Hadi Hashamipour for the help
to provide the figure.

## References

  * (1) G. Aad et al. [ATLAS], “Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC,” Phys. Lett. B 716, 1-29 (2012) [arXiv:1207.7214 [hep-ex]].
  * (2) S. Chatrchyan et al. [CMS], “Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC,” Phys. Lett. B 716, 30-61 (2012) [arXiv:1207.7235 [hep-ex]].
  * (3) S. Chatrchyan et al. [CMS], “Observation of a New Boson with Mass Near 125 GeV in $pp$ Collisions at $\sqrt{s}$ = 7 and 8 TeV,” JHEP 06, 081 (2013) [arXiv:1303.4571 [hep-ex]].
  * (4) M. J. G. Veltman, “The Infrared - Ultraviolet Connection,” Acta Phys. Polon. B 12, 437 (1981).
  * (5) P. Draper and H. Rzehak, “A Review of Higgs Mass Calculations in Supersymmetric Models,” Phys. Rept. 619, 1-24 (2016) [arXiv:1601.01890 [hep-ph]].
  * (6) G. Panico and A. Wulzer, “The Composite Nambu-Goldstone Higgs,” Lect. Notes Phys. 913, pp.1-316 (2016) [arXiv:1506.01961 [hep-ph]].
  * (7) M. Ahmadvand, “Matter and dark matter asymmetry from a composite Higgs model,” Eur. Phys. J. C 81, no.4, 358 (2021) [arXiv:2010.10121 [hep-ph]].
  * (8) B. B. Mandelbrot, “The Fractal Geometry of Nature,” W. H. Freeman and Company: New York, USA (1983).
  * (9) R. Hagedorn, “Statistical thermodynamics of strong interactions at high-energies,” Nuovo Cim. Suppl. 3, 147-186 (1965) CERN-TH-520.
  * (10) S. C. Frautschi, “Statistical bootstrap model of hadrons,” Phys. Rev. D 3, 2821-2834 (1971).
  * (11) G. F. Chew, “Hadron bootstrap hypothesis,” Phys. Rev. D 4, 2330-2335 (1971).
  * (12) I. Bediaga, E. M. F. Curado and J. M. de Miranda, “A Nonextensive thermodynamical equilibrium approach in e+ e- —$>$ hadrons,” Physica A 286, 156-163 (2000) [arXiv:hep-ph/9905255 [hep-ph]].
  * (13) J. Cleymans and D. Worku, “The Tsallis Distribution in Proton-Proton Collisions at $\sqrt{s}$ = 0.9 TeV at the LHC,” J. Phys. G 39, 025006 (2012) [arXiv:1110.5526 [hep-ph]].
  * (14) A. Deppman, “Self-consistency in non-extensive thermodynamics of highly excited hadronic states,” Physica A 391, 6380-6385 (2012) [arXiv:1205.0455 [math-ph]].
  * (15) L. Marques, E. Andrade-II and A. Deppman, “Nonextensivity of hadronic systems,” Phys. Rev. D 87, no.11, 114022 (2013) [arXiv:1210.1725 [hep-ph]].
  * (16) C. Tsallis, “Possible Generalization of Boltzmann-Gibbs Statistics,” J. Statist. Phys. 52, 479-487 (1988).
  * (17) E. M. F. Curado and C. Tsallis, “Generalized statistical mechanics: Connection with thermodynamics,” J. Phys. A 24, L69-L72 (1991) [erratum: J. Phys. A 25, 1019 (1992)] CBPF-NF-041-90.
  * (18) H. Lehmann, K. Symanzik and W. Zimmermann, “On the formulation of quantized field theories,” Nuovo Cim. 1, 205-225 (1955).
  * (19) B. I. Abelev et al. [STAR], “Strange particle production in p+p collisions at s**(1/2) = 200-GeV,” Phys. Rev. C 75, 064901 (2007) [arXiv:nucl-ex/0607033 [nucl-ex]].
  * (20) G. Aad et al. [ATLAS], “Charged-particle multiplicities in pp interactions measured with the ATLAS detector at the LHC,” New J. Phys. 13, 053033 (2011) [arXiv:1012.5104 [hep-ex]].
  * (21) A. Adare et al. [PHENIX], “Identified charged hadron production in $p+p$ collisions at $\sqrt{s}=200$ and 62.4 GeV,” Phys. Rev. C 83, 064903 (2011) [arXiv:1102.0753 [nucl-ex]].
  * (22) K. Aamodt et al. [ALICE], “Production of pions, kaons and protons in $pp$ collisions at $\sqrt{s}=900$ GeV with ALICE at the LHC,” Eur. Phys. J. C 71, 1655 (2011) [arXiv:1101.4110 [hep-ex]].
  * (23) V. Khachatryan et al. [CMS], “Strange Particle Production in $pp$ Collisions at $\sqrt{s}=0.9$ and 7 TeV,” JHEP 05, 064 (2011) [arXiv:1102.4282 [hep-ex]].
  * (24) A. M. Sirunyan et al. [CMS], “Measurement and interpretation of differential cross sections for Higgs boson production at $\sqrt{s}=$ 13 TeV,” Phys. Lett. B 792, 369-396 (2019) [arXiv:1812.06504 [hep-ex]].
  * (25) M. Aaboud et al. [ATLAS], “Combined measurement of differential and total cross sections in the $H\rightarrow\gamma\gamma$ and the $H\rightarrow ZZ^{*}\rightarrow 4\ell$ decay channels at $\sqrt{s}=13$ TeV with the ATLAS detector,” Phys. Lett. B 786, 114-133 (2018) [arXiv:1805.10197 [hep-ex]].
  * (26) F. James, “MINUIT Function Minimization and Error Analysis: Reference Manual Version 94.1,” CERN-D-506.
  * (27) L. Reina, “TASI 2011: lectures on Higgs-Boson Physics,” [arXiv:1208.5504 [hep-ph]].