# Inequality of avalanche sizes in models of fracture

Diksha  Department of Physics, SRM University - AP, Andhra Pradesh 522240,
India  Sumanta Kundu  Department of Physics and Astronomy, University of
Padova, Via Marzolo 8, I-35131 Padova, Italy  INFN, Sezione di Padova, Via
Marzolo 8, I-35131 Padova, Italy  Bikas K. Chakrabarti  Saha Institute of
Nuclear Physics, 1/AF Bidhannagar, Kolkata 700064, India  Soumyajyoti Biswas
Department of Physics, SRM University - AP, Andhra Pradesh 522240, India

###### Abstract

Prediction of an imminent catastrophic event in a driven disordered system is
of paramount importance – from the laboratory scale controlled fracture
experiment to the largest scale of mechanical failure i.e., earthquakes. It
has been long conjectured that the statistical regularities in the energy
emission time series mirrors the “health” of such driven systems and hence
have the potential for predicting imminent catastrophe. Here we use the fiber
bundle model (FBM) and Random fuse model (RFM) of fracture and record the time
series of avalanche sizes. Such avalanches are necessarily very much unequal.
These inequalities in the avalanche time series play an important role in
predicting imminent failure. These inequalities are quantified through some
indices, namely the Hirsch-index $h$, Gini index $g$, and Kolkata index $k$
that are estimated from the Lorenz function $L(p)$ where the $p$ fraction of
avalanche size acquires the cumulative fraction of avalanche mass when
arranged in the ascending order of their sizes. Here we calculate these
inequalities indices using the time series of avalanche size. In the mean-
field case of the FBM, we can calculate $h,g,$ and $k$ analytically: these
inequalities’ endpoints or terminal values can be a very good indicator of
imminent catastrophic failure. For non-mean field cases, we give numerical
results. These inequalities of avalanche size show universal scaling
properties near the failure point.

##  I Introduction

When a solid is slowly loaded until its breaking point, it goes through many
microscopic damages that grow, interact and accumulate into a macroscopic
failure (for example, a shear plane in the case of compression induced damage
in quasi-brittle objects). However, all breaking events are not equally
damaging. Indeed, if we go by the acoustic energy signal ($S$) emitted in each
event, their size distribution ($P(S)$) usually has a fat-tail ($P(S)\sim
S^{-\delta}$), implying a very few events accounting for most of the damages
incurred Biswas  _et al._ (2015); Bonamy and Bouchaud (2011). This ‘damage
inequality’ is reminiscent of inequalities observed in a myriad of complex
systems, including physical (from crackling noise Sethna  _et al._ (2001) in
magnetic domain wall movements to the largest mechanical failures in
earthquakes Gutenberg and F. Richter (1949)), sociological Chatterjee  _et
al._ (2015)(deaths in armed conflicts, citations of papers) and economic
(wealth distribution) Chakrabarti  _et al._ (2013) systems. This fat-tail or
power-law behavior is often attributed to the vicinity of a critical point in
the system that is reached either through the fine tuning of an external
parameter (magnetic fields in the case of domain walls) or the system self-
organizes to the critical point through a slow drive Bak  _et al._ (1987).

It was empirically noted by Pareto in 1896 in his famous 80-20 law that 80% of
wealth is accumulated with 20% of people. Since then this observation was much
broadly applied to see that in many socio-economic systems 80% of ‘successes’
come from 20% of ‘attempts’ Pareto (1896-1965). Since then, much has been
achieved in studies of inequality, primarily in socio-economic systems, mostly
because of the immense adverse effect such inequalities could have in those
contexts Bourguignon (2015). A much recent interest is to understand the
inequality in the responses of physical systems, particularly when such
systems are near a critical point Biswas and Chakrabarti (2021). The main goal
in such a scenario is to forecast imminent catastrophe (or distance to the
critical point) in such systems. For example, it was shown that in models of
fracture, the growing inequality of avalanche sizes shows universal behavior
near critical point Biswas and Chakrabarti (2021). It was then shown that in a
supervised machine learning algorithm, measures of inequality were the most
important attribute in predicting the imminent failure point Diksha and Biswas
(2022). It was also shown subsequently that a family of self-organized
critical (SOC) systems (e.g., BTW and Manna sandpile models among others) also
show universal behavior in terms of the inequality measures of their responses
(avalanches) Manna  _et al._ (2022). Along that line, it was shown using
extensive data analysis in a wide variety of socio-economic systems that have
been conjectured to be in self-organized critical states, that the measures of
inequalities in such systems show the similar behavior as seen in the SOC
models. Such observations underline the importance of the inequality measures
of responses of a system near its critical point for uncovering the formal
link between such measures and critical phenomena Ghosh  _et al._ (2022) and
also to make use of such a link for the practical purposes of forecasting an
imminent critical point in systems where such a point could have a
catastrophic consequence (for example, in fractures and earthquakes).

In this work, we study two generic models of fracture in disordered solids –
the fiber bundle model (FBM) and the random fuse model (RFM) in terms of the
inequality of their response functions (avalanche time series) as they are
stressed to their breaking points. Particularly, in the mean-field FBM, we
analytically calculate the inequality indices and their scaling behavior near
the critical point, thereby illustrating the previously numerically observed
universal behavior. We also study the more realistic non-mean field limit of
the FBM and the RFM and show numerically that similar scaling are observed
(with different exponent values) in such cases as well.

The paper is organized as follows: first, we define the indices that were used
to quantify the inequalities in the avalanche statistics of these models, then
we describe the models and their simulations. We then go on to present the
analytical calculations for the mean-field case and compare them with previous
numerical studies. Then we give the numerical results for the non-mean-field
case of the FBM and the RFM. Finally, we discuss the results and their
implications for forecasting imminent failure in these models and potentially
in experiments and conclude.

##  II Quantification of inequality of avalanches: The inequality indices

The inequalities in socio-economic contexts have been quantified using several
indices, for example, the Gini-index $g$ Dalton (1920) (for wealth), the
Hirsch-index $h$ (for citations) Hirsch (2005), and more recently the Kolkata
index $k$ Ghosh  _et al._ (2014), to name a few. The role of these indices in
parametrizing socio-economic inequality has expanded over the past few
decades. These inequality measures are defined through the Lorenz curve. When
a time series is arranged in ascending order of the size of the values in the
series, the Lorenz function $L(p)$ gives the cumulative fraction of the total
‘mass’ acquired by the $p$ fraction of the smallest events. If all the values
are equal in size, then the Lorenz function would be a diagonal line from the
origin to (1,1), known as the equality line. But in general, these values are
not equal in size, so the Lorenz function is nonlinear, always staying under
the equality line and increasing monotonically from $L(p=0)=0$ to $L(p=1)=1$.
The Gini index is then given by twice the area between the Lorenz curve and
the equality line, where $g=0$ implies perfect equality and $g=1$ shows the
extreme inequality. The crossing point of the opposite diagonal line with the
Lorenz function gives the value of the Kolkata index $k$, which says that
$(1-k)$ fraction of the largest events account for the $k$ fraction of the
total ‘mass’, where $k=0.5$ shows the perfect equality and $k=1$ shows the
extreme inequality. The Hirsch-index is generally calculated in the case of
citations of individual scientists or scholars but the definition could be
expanded to measure the $h$ index of any series. The $h$-index is the highest
number $h$ such that $h$ events have at least a size $h$ each. For a more
detailed discussion, see ref. Ghosh  _et al._ (2022)

In this work, we will be measuring the inequality of the time series of the
avalanches recorded for slowly loaded disordered materials, represented by
generic threshold activated models. Therefore, for our purposes, the Lorenz
function $L(p)$ would represent the fraction of the breaking events resulting
from the smallest $p$ fraction of the avalanches. The other definitions would
follow the same line. We will show how the inequalities grow as the system
gradually approach the catastrophic failure points. Most importantly, we will
monitor the values ($g_{f},k_{f},h_{f}$) the inequality indices take just
prior to the catastrophic breakdown. We will be interested in the universality
of these values and also the finite size and off-critical scaling behavior
these measures show, particularly for the context relevant to forecasting as
an imminent failure.

##  III Models and Methods

As mentioned before, in this work we consider two generic models of fracture –
the fiber bundle model (FBM) T. Pierce (1926); Pradhan  _et al._ (2010) and
the random fuse model (RFM) de Arcangelis, L.  _et al._ (1985); Kahng  _et
al._ (1988); Zapperi  _et al._ (1997). In the FBM, an ensemble of $N$ linear
elastic fibers are considered, which have different individual breaking
thresholds. When a load is applied on all the fibers, some of the weak fibers
fail and the load carried by those fibers are redistributed among the
remaining ones, which might in-turn fail, leading to an avalanche event. At a
critical value of the load, the entire system collapses. In the RFM, an
ensemble of electrical fuses are arranged in a network (say, a square
lattice). Each fuse has a failure threshold, beyond which it is burnt. When a
voltage difference is applied across the system (say, the two ends of a square
lattice), the weaker fuses are burnt, requiring the current to be diverted
through other paths, which might in-turn cause other fuse burning, leading to
an avalanche event. The modification in the flow of current due to burnt fuses
mimics the modification of the stress field in an elastic solid due to local
damages. Both of these models are very well studied in the literature from
various viewpoints of applications in fracture of disordered systems Alava
_et al._ (2006); Hansen  _et al._ (2015).

Finding an indication of the imminent failure in these models, without using
detailed information regarding the individual failure thresholds, therefore,
is a long-standing crucial issue. This has been approached through various
different ways, including changes in the avalanche size distribution exponent
prior to failure Pradhan  _et al._ (2005), non-monotonic behavior of the
elastic energy stored in the system Pradhan  _et al._ (2019) and so on.
However, here we will study the behavior of the inequality indices, which have
shown promising indications as precursory measures of imminent breakdown
Diksha and Biswas (2022).

###  III.1 Simulation methods

While a mean-field limit of the FBM is analytically tractable, the more
realistic versions of FBM and also the RFM are not. Those limits are accessed
through numerical simulations. Therefore, here we describe the simulation
methods for the two models using which the avalanche time series could be
obtained. The inequality indices are then measured on those time series.

####  III.1.1 Simulating FBM

In the FBM, a large number of parallel fibers are connected between two
horizontal plates. The top plate with hanging fiber is rigid and the rigidity
of the lower plate determines the interaction range of the model. All the
fibers are here assumed to have identical elastic constants but have a
different failure threshold from each other. Here the failure thresholds of
the fibers are randomly assigned from the uniform distribution between (0:1).
When a load is applied to the system, the weakest fibers fail irreversibly and
the load of the broken fibers is redistributed to the remaining intact fibers.
The increased load on the surviving fibers then may cause more failures and
this process will go on until the remaining fiber are strong enough to hold
the extra load or the whole bundle collapses. The external load is kept
constant throughout the redistribution process and it is again increased to
initiate the dynamics once the dynamics stops. the number of fibers breaking
between two successive load increments is called an avalanche. The
redistribution of the load mentioned above depends on the elastic properties
of the lower plate. The two cases of load redistribution are equal load
sharing (ELS) and Local load sharing (LLS). In ELS, when the lower plate is
absolutely rigid, the load of the broken fiber is equally redistributed to all
surviving fibers. In LLS, when the lower plate deforms under loading, then a
higher fraction of the load of broken fiber is carried by the surviving
neighbors that are nearer to the broken fiber. The extreme case is the nearest
neighbor load sharing. But a long-range (power-law) load sharing has also been
studied before Hidalgo  _et al._ (2002). Here we redistribute the load of the
broken fiber between $R$ surviving neighbors on either side (see e.g., Roy
_et al._ (2017)), when the fibers are arranged in an one-dimensional lattice
of length $L$. Clearly, $R=1$ is the nearest neighbor load sharing, and $R\sim
L$ is the ELS limit (in fact the ELS limit is reached much earlier Roy  _et
al._ (2019)).

Here we study both the mean-field limit (ELS) and the local load sharing limit
with finite values of $R$. We mention the respective cases in the results as
appropriate.

####  III.1.2 Simulating RFM

The random fuse network model that we consider here is defined on a two-
dimensional tilted square lattice of linear size $L$ with periodic boundary
conditions along the horizontal direction. A potential difference is applied
along the vertical direction between the two opposite sides of the lattice.
Each bond $i$ representing an Ohmic resistor with unit conductance carries
electrical current $i_{i}$ until it burns out irreversibly at a threshold
value $b_{i}$ of its current and then, it becomes an insulator. The thresholds
values $\\{b_{i}\\}$ are selected randomly and independently from a power-law
distribution: $p(b)\sim b^{-1}$, bounded between $10^{-\beta}$ and
$10^{\beta}$. This is a generic form of distribution with a decaying power-law
tail that has been considered widely for modeling fracture of heterogeneous
materials with varying degrees of disorder Moreira  _et al._ (2012); Roy  _et
al._ (2015); Hope  _et al._ (2015); Roy  _et al._ (2017).

Initially, all the bonds are intact and the applied voltage difference $V$ is
raised quasi-statically from zero. The specific geometry of our lattice
ensures that every bond initially carries the same amount of current.
Therefore, the breaking process initiates when the current through the
“weakest bond” in the system reaches its breaking threshold at
$V=V_{0}=(L-1)\times\min\\{b_{i}\\}$. Subsequently, the bond burns and it is
irreversibly removed. The new current distribution is then determined by
keeping the external voltage difference fixed at $V=V_{0}$. This may initiate
a sequence of bond burning. At this stage, all the bonds carrying currents
higher than their respective threshold values are removed simultaneously and
again, the current distribution is calculated. Numerically, a set of
Kirchhoff’s equations are solved using the conjugate gradient method to
determine the current distribution George Batrouni and Hansen (1988). This way
the breaking process continues until a stable state is reached when the
current through all the remaining intact bonds is lower than their respective
breaking thresholds. This completes an avalanche. The largest ratio between
current and threshold, i.e., $\max(i_{i}/b_{i})$ is then calculated for all
intact bonds $i$ to determine the next weakest bond to be burnt and removed.
This initiates a new avalanche. Accordingly, the external voltage difference
is raised.

In this voltage-controlled setting, the breaking process stops completely when
a final crack comprising of burnt bonds wraps around the lattice in a
direction transverse to the applied potential difference, i.e., no current
flows through the system. The number of bonds that burn between two successive
increments of the external voltage determines the avalanche size.

Note that at an early stage of the breaking process, the intact bond having
the minimum value of the breaking threshold determines the weakest bond.
However, depending on the value of $\beta$, this is not generally true at the
later stage. The current concentration around the crack zone competes with the
local strength of the breaking thresholds. As a consequence, the breaking
events become correlated Zapperi  _et al._ (1997); Moreira  _et al._ (2012);
Shekhawat  _et al._ (2013); Hansen and Schmittbuhl (2003). Certainly, the
choice of $\beta$ is important as discussed in Ref. Moreira  _et al._ (2012).
In the limit of weak disorder, i.e., $\beta\to 0$, the width of the
distribution is so narrow that the fracture process becomes localized and only
a single crack grows in the system, while it becomes percolation-like in the
limit of strong disorder, i.e., $\beta\to\infty$. In the later case, one only
gets avalanches of size unity and always burns the intact bond with the
smallest breaking threshold on the conducting backbone. Importantly, depending
on $L$, the breaking dynamics is avalanche dominated in the intermediate range
of disorder Moreira  _et al._ (2012); Shekhawat  _et al._ (2013).

We choose the value of $\beta=$1 and consider system sizes $L=$ 32, 64, 128,
256, and 512 such that the system remains in the transition regime from the
weak to the strong disorder where avalanche dynamics can be observed.

##  IV Results

We now move on to describe the behavior of the inequality indices $g$, $h$,
and $k$ for the avalanche statistics of the FBM and the RFM. We first present
the analytical results for the FBM in the mean-field (ELS) limit and then go
on to discuss the local load sharing version and the RFM results.

###  IV.1 Inequality indices for FBM in mean field limit: Analytical results

The FBM in the mean-field (ELS) limit is analytically tractable Pradhan  _et
al._ (2010). The system size (initial number of fibers) is denoted by $N$, the
applied load is $W$, giving the load per fiber $\sigma=W/N$, with $\sigma_{c}$
denoting the critical value at which catastrophic breakdown happens. The
fraction of intact fibers at load per fiber $\sigma$ is given by

$U(\sigma)=U^{*}(\sigma_{c})+(\sigma_{c}-\sigma)^{1/2}.$ (1)

Figure 1: The plot of $k(q)$ (obtained numerically from Eq. (16)) is against
$g(q)$ (obtained from Eq. (15)) for the mean field FBM. The line
$k=\frac{1}{2}+\frac{3}{8}g$ is also plotted for reference.

The magnitude of the change in this fraction for a differential increment in
the load per fiber value is, therefore,

$S(\sigma)=\left|\frac{dU(\sigma)}{d\sigma}\right|=\frac{1}{2}(\sigma_{c}-\sigma)^{-1/2}.$
(2)

This is a measure of avalanche size at $\sigma$, when the load increment is
done by a fixed amount $\delta=d\sigma$ at every step.

####  IV.1.1 Terminal values of the inequality indices

The function in Eq. (2) is monotonically increasing up to the failure point
$\sigma=\sigma_{c}$, hence the Lorenz function at the point of collapse of the
system is given by

$L_{f}(p)=\frac{\int\limits_{0}^{p\sigma_{c}}(\sigma_{c}-\sigma)^{-1/2}d\sigma}{\int\limits_{0}^{\sigma_{c}}(\sigma_{c}-\sigma)^{-1/2}d\sigma}$
(3)

Figure 2: Finite-size scaling of $g$ for system sizes ($N/10^{3}$ = 5, 30, 50,
60, 80, 100, 200, 300, 400, 500, 700) for the mean-field fiber bundle model.
From this figure, we can see that $g_{f}$ shows the finite size scaling for
the larger system size.

On simplification, this gives

$L_{f}(p)=1-\sqrt{1-p},$ (4)

where $\sigma_{c}$ cancels out, giving the function the robustness it shows
against various threshold distributions (see Fig. (11) in ref. Diksha and
Biswas (2022)). The Gini index ($g$) at the point of collapse is then

$g_{f}=1-2\int\limits_{0}^{1}L_{f}(p)dp=\frac{1}{3}.$ (5)

The Kolkata index ($k$) at the point of collapse can be found from

$1-k_{f}=L_{f}(k_{f})=1-\sqrt{1-k_{f}},$ (6)

which gives

$k_{f}=\frac{\sqrt{5}-1}{2},$ (7)

which is the inverse golden ratio. The other solution is irrelevant.

It is also known Pradhan  _et al._ (2010) that for the fiber bundle model with
an equal amount of load increase in each step, the number of fibers broken due
to the increase of load between $m\delta$ to $(m+1)\delta$ is given by

$S(m)=\frac{\delta}{\sqrt{1-4m\delta/N}},$ (8)

where $m=0,1,2,\dots,N/4\delta$. This suggests a monotonic increase in the
avalanche size, which is strictly true for equispaced thresholds (given by
$i/N$, $i=1,\dots,N$) but the scaling relations are valid for non-equispaced
thresholds as well. To find the Hirsch index ($h$), the avalanches should be
first sorted in descending order, which we do with the transformation
$m\to\frac{N}{4\delta}-m$, giving

$S\left(\frac{N}{4\delta}-m\right)=\frac{\delta}{\sqrt{1-\frac{4\delta}{N}\left(\frac{N}{4\delta}-m\right)}}=\frac{\delta}{\sqrt{\frac{4\delta
m}{N}}}.$ (9)

Then the value of the Hirsch index at the point of collapse can be found from

$\frac{\delta}{\sqrt{\frac{4\delta h_{f}}{N}}}=h_{f},$ (10)

which gives

$h_{f}\propto N^{1/3}.$ (11)

This scaling is in contrast with the previously conjectured form (see ref.
Biswas and Chakrabarti (2021)). However, as we shall shortly see, this is
obeyed in accurate numerical simulations of the model.

####  IV.1.2 Off-critical scaling

Now, if the loading is done upto $\sigma=q\sigma_{c}$, where $q\leq 1$, then
the Lorenz function can be found from

$L(p,q)=\frac{\int\limits_{0}^{pq\sigma_{c}}(\sigma_{c}-\sigma)^{-1/2}d\sigma}{\int\limits_{0}^{q\sigma_{c}}(\sigma_{c}-\sigma)^{-1/2}d\sigma},$
(12)

which gives

$L(p,q)=\frac{1-\sqrt{1-pq}}{1-\sqrt{1-q}},$ (13)

which, in the limit $q=1$, gives back $L_{f}$ in Eq. (4). It is then possible
to calculate the evolution of $g$ and $k$ as the loading is increased, using

$g(q)=1-2\int\limits_{0}^{1}L(p,q)dp,$ (14)

which gives on simplification

$g=1-\frac{4(1-q)\sqrt{1-q}+6q-4}{3q(1-\sqrt{1-q})}.$ (15)

Similarly, the equation for $k(q)$ reads

$k(q)=1-\frac{1-\sqrt{1-k(q)q}}{1-\sqrt{1-q}}.$ (16)

For any value of $q$, the value of $k(q)$ can be numerically estimated with
arbitrary accuracy from the equation above. In Fig. 1, the values obtained for
$k(q)$ from the above equation is plotted against those obtained for $g(q)$
from Eq. (16). A phenomenologically derived expression
$k=\frac{1}{2}+\frac{3}{8}g$ Ghosh  _et al._ (2022) is also plotted for
reference. It shows good agreement for lower values of $g$, for which the
phenomenological argument is valid.

Figure 3: The finite size scaling of h-index, following Eqs. (22) & (23) for
the mean field FBM. Figure 4: The plot of $g(q)$ against $k(q)$ for different
interaction ranges $(R=1,100,500,1000,3000,5000,7000)$ for three different
system sizes $(L=10000,30000,50000)$ in LLS fiber bundle model. The plot shows
that terminal values of $g$ and $k$ are independent of the system sizes. With
increasing the value of $R$, the terminal values of $g$ and $k$ approach to
the values obtained analytically in the mean-field limit.

From Eq.(15), the $q\to 1$ limit ($\sigma\to\sigma_{c}$) yields

$g-g_{f}\propto\sqrt{1-q},$ (17)

meaning that we could write all response functions in terms of the interval
$|g-g_{f}|$ instead of $\sigma-\sigma_{c}$ with the appropriate changes in the
scaling exponent (see ref.Das and Biswas (2023) for details). An important
consequence of the above equation is the finite size scaling of the terminal
value of $g$. We could define a correlation length exponent as
$\xi\propto|\sigma-\sigma_{c}|^{-\nu}$ Biswas and Chakrabarti (2020), which
would span the system size at the critical point. We should then have, for a
finite size system at its critical point,
$|\sigma_{c}(N)-\sigma_{c}(N\to\infty)|\propto N^{-1/\nu}$. Also, at the
critical point of the finite system, we should have from Eq. (17),

$g_{f}(N)-g_{f}(N\to\infty)\propto|\sigma_{c}(N)-\sigma_{c}(N\to\infty)|^{1/2}\propto
N^{-1/2\nu}.$ (18)

Given that $\nu=3/2$ Roy  _et al._ (2015); Biswas and Chakrabarti (2020), we
should have $|g_{f}(N)-g_{f}(N\to\infty)|\propto N^{-1/3}$, which is what we
see in Fig. 2.

To find the load and size dependence of $h$, consider again the Eq. (8).
Suppose we stop the loading at $m=m_{max}=\frac{qN}{4\delta}$, where $q\leq
1$, as before. The avalanches arranged in descending order would then read

$S(m_{max}-m)=\frac{\delta}{\sqrt{1-\frac{4\delta}{N}(m_{max}-m)}}=\frac{\delta}{\sqrt{1-\frac{4\delta}{N}(\frac{qN}{4\delta}-m)}}.$
(19)

The Hirsch index can then be found, as before, from

$\frac{\delta}{\sqrt{1-\frac{4\delta}{N}(\frac{qN}{4\delta}-h)}}=h.$ (20)

This can be simplified to

$h^{3}+\frac{N(1-q)}{4\delta}h^{2}-\frac{N\delta}{4}=0.$ (21)

This cubic equation can be solved. But the expression for $h$ is then too
complex to figure out the scaling form near the critical point. Instead, we
can approximate that near the critical point ($q\to 1$), the cubic term in $h$
would be the most dominant, hence responsible for the leading order scaling.
In that case, the quadratic term can be replaced by its value exactly at the
critical point $h_{f}^{2}=\left(\frac{\delta N}{4}\right)^{2/3}$. Then, near
the critical point $\sigma_{c}=1/4$,

$\frac{4\delta}{N}h^{3}\approx\delta^{2}-(1-q)\left(\frac{\delta
N}{4}\right)^{2/3}.$ (22)

This gives

$h^{3}\approx\frac{N\delta}{4}-\frac{1-q}{4}\frac{N}{\delta}\left(\frac{\delta
N}{4}\right)^{2/3}.$ (23)

On simplification, this gives

$\displaystyle h^{3}$ $\displaystyle\approx$ $\displaystyle\frac{\delta
N}{4}-\frac{1-q}{4}\frac{N}{\delta}\left(\frac{\delta N}{4}\right)^{2/3}$ (24)
$\displaystyle\approx$ $\displaystyle
N\left[\frac{\delta}{4}-(\sigma_{c}-\sigma)N^{2/3}\left(\frac{1}{4\sqrt{\delta}}\right)^{2/3}\right],$

where $\frac{1-q}{4}=\sigma_{c}-\sigma$ (remembering $\sigma=q\sigma_{c}$).
This is now in the scaling form

$h\approx N^{1/3}f\left[(\sigma_{c}-\sigma)N^{1/\nu}\right],$ (25)

with $\nu=3/2$ Roy  _et al._ (2015) being the correlation length exponent.

Fig. 3 shows the numerical validation of the scaling form above. The
simulation was done for a fixed load increase with $\delta=10$. So, as
$\sigma\to\sigma_{c}$, the quantity $h/N^{1/3}$ should tend to
$(5/2)^{1/3}\approx 1.357$, which it does (as can be seen from the inset of
Fig. 3). Also, away from the critical point, the quantity should be
proportional to $(\sigma_{c}-\sigma)^{-1/2}$, which is also seen. The data
collapse for the three orders of magnitude system size variation confirms the
scaling $h_{f}\sim N^{1/3}$.

###  IV.2 Inequality indices in the non-mean field limit (local load sharing
case) of FBM: Simulation results

This is the case when a load of a broken fiber is redistributed to the nearest
$R$ surviving fibers when they are arranged in a one-dimensional lattice of
size $L$. Here we take the threshold distribution to be uniform between 0 and
1. Here $R=1$ means, the sharing of the load is fully local. In this case, we
take three different system sizes ($L=10000,30000,50000$) and varied the value
of $R$, and then calculated the value of $h$, $g$, and $k$. In Fig. 4 we plot
the $g$ against $k$ and we have observed numerically when the interaction
range ($R$) is increased, $g_{f}$ and $k_{f}$ values are very close to the
analytical results obtained for the mean-field case. Prior to reaching the
terminal values, the variation of $g$ vs $k$ almost follows the mean field
variation, except that the growth stops earlier (usually there is less number
of avalanches for local load sharing). In Fig. 5 the variation of $h$ is
plotted for different values of $R$. Here also, the terminal values are
smaller for lower values of $R$, and so are their system size scaling. Indeed,
from the Fig. 6 it can be seen that for fully local load sharing ($R=1$),
scaling of $h_{f}$ shows the logarithmic dependence with $L$, but as we
increase the interaction range, the scaling of $h_{f}$ shows the power law
dependence of $L^{1/3}$, similar to what we expect in the mean-field limit.

Figure 5: The variation of h-index with re-scaled time $t/L$ for different
interaction ranges $(R=1,5,10,500,1000,3000)$ for three different system sizes
$(L=10000,30000,50000)$ for LLS fiber bundle model. Figure 6: System size $L$
dependence of the terminal values of $h=h_{f}$ for different interaction
ranges $R$ for one dimensional LLS fiber bundle model in a log-log scale. It
is seen that for $R=1$ when the load sharing is fully local $h_{f}$ shows the
logarithmic dependence, and when $R=$ 100, and 500, the $h_{f}$ shows the
power law dependence. Here the solid blue and black lines indicate the
$\log(L)$ and $L^{1/3}$ functional form, respectively.

###  IV.3 Inequality indices for the RFM: Simulation results

Finally, we compare the above results for a more realistic model of fracture,
the random fuse model (RFM). In this case, the threshold distribution is a
power law with exponent $-1$, as mentioned before. The limits of the
distribution are between $10^{-\beta}$ to $10^{\beta}$. We choose the value of
$\beta$=1 for the system sizes $L=$ 32, 64, 128, 256 & 512\. We calculated the
inequality measures using 500 samples for system sizes up to $L=256$ and 40
samples for $L=512$. Given that $N=L\times L$, if we expect a similar scaling
as in the case of FBM, we should have $h_{f}\sim L^{2/3}$. In Fig. 7, the
function $h/L^{2/3}$ collapses on a single curve, indicating a similar scaling
as that in the FBM.

The values of $g_{f}$ and $k_{f}$ are set by the avalanche size distribution
exponent. Particularly, if the avalanche size distribution exponent is
$\delta$, then it is straightforward to see that when the avalanches are
arranged in the ascending order of their sizes, the resulting curve will
diverge with an exponent $n=1/(\delta-1)$. For $\beta=1.0$, $\delta\approx
1.99$. Now, it is also possible to show Das and Biswas (2023) that for $n>1$,
$g_{f}=k_{f}=1$. Therefore, in the finite size scaling of their values, we
would have to put $g_{f}(L\to\infty)\to 1$, resulting in a scaling form as
before $|1-g_{f}(L)|\propto L^{-1/2\nu}$. In Fig. 8, we do see a power law
decay of $|1-g_{f}(L)|$ with $L$, but not with the expected exponent value
(typically $\nu=4/3$ for RFM Ramstad  _et al._ (2004)). Nevertheless, the
power law scaling is seen, with an exponent $-0.23$.

Figure 7: The $h$-index (scaled by $L^{2/3}$) is plotted against the time
(scaled by time/L) for different system sizes $(L=32,64,128,256)$ for the two-
dimensional random fuse model with power-law distributed breaking threshold.
The inset shows the variation of $h$ with time for the different system sizes.

##  V Discussion and Conclusion

Near a critical point, the responses of a system show critical divergence.
This necessarily means that the responses are highly unequal as a critical
point is gradually approached. In the case of disordered solids and their
models, these unequal responses are characterized in terms of the avalanche
time series, which are accumulated local breaking events, experimentally
detected as acoustic emissions. It is well known that the statistical
properties of such acoustic emissions are useful to characterize the ‘health’
of a system i.e., to estimate the vicinity of an approaching critical
(breaking) point.

Figure 8: The $g$-index is plotted against time for different system sizes
$(L=32,64,128,256,512)$ for the random fuse model with power-law distributed
breaking threshold. As can be seen from the figure the terminal values of g is
changing with the system size. The inset shows the variation of the terminal
values of the g-index with time in the log-log scale, exhibiting a power-law
scaling.

Recently, quantitative measures of the inequality of such avalanches, e.g.,
the Hirsch index ($h$), the Gini index ($g$), and the Kolkata index ($k$) have
been used to predict the imminent critical point in physical systems Biswas
and Chakrabarti (2021). Particularly, it was noted that $g$ and $k$ reach
universal terminal values $g_{f}$ and $k_{f}$ at the critical point,
irrespective of the details of the system (types of disorder present in them).
It is therefore very useful to monitor the inequality of the avalanche
statistics to estimate the distance from an approaching critical point Diksha
and Biswas (2022). It is more useful to monitor the inequality indices than
the other universal quantities e.g., the avalanche size distribution critical
exponent because such exponent shows universal properties only in the
asymptotic limit of the large avalanche sizes, where the failure is already
extremely imminent. Indeed, extensive numerical simulations and machine
learning studies have indicated the same Diksha and Biswas (2022).

Here we have shown analytically, in the mean-field FBM, the scaling behavior
of the critical indices ($h$, $k$, and $g$). Particularly, for the terminal
values of the inequality indices we have shown that $h_{f}\sim N^{1/3}$, where
$N$ is the system size, $g_{f}=1/3$ and $k_{f}=(\sqrt{5}-1)/2$. We have also
shown the off-critical Widom-Stauffer like scaling for $h$ (see Eq (23)).
These are also verified numerically (see Figs 3 & 6). Some of these scaling
relations differ from what was conjectured using only numerical simulations
before Biswas and Chakrabarti (2021).

We then study the more realistic local load sharing non-mean field version of
the FBM and also the RFM. In LLS-FBM, we see that the scaling of $h_{f}$
shifts from a logarithmic dependence for very localized load sharing, to the
power-law dependence mentioned above, as the load sharing range $R$ is
increased. However, the power-law scaling with system size is observed for the
RFM (see Fig 7). We also show the finite size effect in $g_{f}$ and $k_{f}$
for the FBM and the RFM (see Figs 2 & 8 ) and note that for the Local-Load-
Sharing FBM and RFM models having avalanche size exponent $\delta<2$, the
terminal values ($g_{f}$) of $g$ and ($k_{f}$) of $k$ both approach unity as
the system size $L$ approaches infinity.

In conclusion, the inequalities of the avalanche sizes in models of fracture
show universal scaling properties near the failure point. Here we show
analytically (for mean field FBM) and numerically (for non-mean field FBM and
RFM) the scaling properties of the inequality indices. These results show that
monitoring the inequality indices in stressed disordered solids can indicate
approaching a critical point, as was conjectured from numerical and machine
learning analysis before.

## acknowledgement

SK is supported by the research grant ‘ORLA$\\_$BIRD2020$\\_$01’ of the
University of Padova. The simulations of FBM were done using HPCC Surya in SRM
University - AP.

## References

  * Biswas  _et al._ (2015) S. Biswas, P. Ray, and B. K. Chakrabarti,  _Statistical physics of fracture, breakdown and earthquake_ (Wiley, 2015).
  * Bonamy and Bouchaud (2011) D. Bonamy and E. Bouchaud, “Failure of heterogeneous materials: A dynamic phase transition?” Physics Reports 498, 1–44 (2011).
  * Sethna  _et al._ (2001) James P. Sethna, Karin A. Dahmen, and Christopher R. Myers, “Crackling noise,” Nature 410, 242–250 (2001).
  * Gutenberg and F. Richter (1949) B. Gutenberg and C. F. Richter, “Seismicity of the earth and associated phenomena,” Princeton University Press , 295 (1949).
  * Chatterjee  _et al._ (2015) Arnab Chatterjee, Asim Ghosh, Jun ichi Inoue, and Bikas K Chakrabarti, “Social inequality: from data to statistical physics modeling,” Journal of Physics: Conference Series 638, 012014 (2015).
  * Chakrabarti  _et al._ (2013) B.K. Chakrabarti, A. Chakraborti, S.R. Chakravarty, and A. Chatterjee,  _Econophysics of Income and Wealth Distributions_ (Cambridge University Press, 2013).
  * Bak  _et al._ (1987) Per Bak, Chao Tang, and Kurt Wiesenfeld, “Self-organized criticality: An explanation of the 1/f noise,” Phys. Rev. Lett. 59, 381–384 (1987).
  * Pareto (1896-1965) V. Pareto, “Cours d’ economie politique,” Reprinted as a volume of Oeuvres Compl‘etes (1896-1965).
  * Bourguignon (2015) F. Bourguignon, “Globalization of inequality,” Princeton University Press (2015).
  * Biswas and Chakrabarti (2021) Soumyajyoti Biswas and Bikas K. Chakrabarti, “Social inequality analysis of fiber bundle model statistics and prediction of materials failure,” Phys. Rev. E 104, 044308 (2021).
  * Diksha and Biswas (2022) Diksha and Soumyajyoti Biswas, “Prediction of imminent failure using supervised learning in a fiber bundle model,” Phys. Rev. E 106, 025003 (2022).
  * Manna  _et al._ (2022) S.S. Manna, Soumyajyoti Biswas, and Bikas K. Chakrabarti, “Near universal values of social inequality indices in self-organized critical models,” Physica A: Statistical Mechanics and its Applications 596, 127121 (2022).
  * Ghosh  _et al._ (2022) Asim Ghosh, Soumyajyoti Biswas, and Bikas K. Chakrabarti, “Success of social inequality measures in predicting critical or failure points in some models of physical systems,” Frontiers in Physics 10 (2022), 10.3389/fphy.2022.990278.
  * Dalton (1920) Hugh Dalton, “The measurement of the inequality of incomes,” The Economic Journal 30, 348–361 (1920).
  * Hirsch (2005) J. E. Hirsch, “An index to quantify an individual’s scientific research output,” Proceedings of the National Academy of Sciences 102, 16569–16572 (2005).
  * Ghosh  _et al._ (2014) Asim Ghosh, Nachiketa Chattopadhyay, and Bikas K. Chakrabarti, “Inequality in societies, academic institutions and science journals: Gini and k-indices,” Physica A: Statistical Mechanics and its Applications 410, 30–34 (2014).
  * T. Pierce (1926) F. T. Pierce, “32—x.—tensile tests for cotton yarns v.—“the weakest link” theorems on the strength of long and of composite specimens,” Journal of the Textile Institute Transactions 17, T355–T368 (1926).
  * Pradhan  _et al._ (2010) Srutarshi Pradhan, Alex Hansen, and Bikas K. Chakrabarti, “Failure processes in elastic fiber bundles,” Rev. Mod. Phys. 82, 499–555 (2010).
  * de Arcangelis, L.  _et al._ (1985) de Arcangelis, L., Redner, S., and Herrmann, H. J., “A random fuse model for breaking processes,” J. Physique Lett. 46, 585–590 (1985).
  * Kahng  _et al._ (1988) B. Kahng, G. G. Batrouni, S. Redner, L. de Arcangelis, and H. J. Herrmann, “Electrical breakdown in a fuse network with random, continuously distributed breaking strengths,” Phys. Rev. B 37, 7625–7637 (1988).
  * Zapperi  _et al._ (1997) Stefano Zapperi, Purusattam Ray, H. Eugene Stanley, and Alessandro Vespignani, “First-order transition in the breakdown of disordered media,” Phys. Rev. Lett. 78, 1408–1411 (1997).
  * Alava  _et al._ (2006) Mikko J. Alava, Phani K. V. V. Nukala, and Stefano Zapperi, “Statistical models of fracture,” Advances in Physics 55, 349–476 (2006), https://doi.org/10.1080/00018730300741518 .
  * Hansen  _et al._ (2015) Alex Hansen, Per Christian Hemmer, and Srutarshi Pradhan,  _The Fiber Bundle Model: Modeling Failure in Materials_ (Wiley, 2015).
  * Pradhan  _et al._ (2005) Srutarshi Pradhan, Alex Hansen, and Per C. Hemmer, “Crossover behavior in burst avalanches: Signature of imminent failure,” Phys. Rev. Lett. 95, 125501 (2005).
  * Pradhan  _et al._ (2019) Srutarshi Pradhan, Jonas T. Kjellstadli, and Alex Hansen, “Variation of elastic energy shows reliable signal of upcoming catastrophic failure,” Frontiers in Physics 7 (2019), 10.3389/fphy.2019.00106.
  * Hidalgo  _et al._ (2002) Raul Cruz Hidalgo, Yamir Moreno, Ferenc Kun, and Hans J. Herrmann, “Fracture model with variable range of interaction,” Phys. Rev. E 65, 046148 (2002).
  * Roy  _et al._ (2017) Subhadeep Roy, Soumyajyoti Biswas, and Purusattam Ray, “Modes of failure in disordered solids,” Phys. Rev. E 96, 063003 (2017).
  * Roy  _et al._ (2019) Subhadeep Roy, Soumyajyoti Biswas, and Purusattam Ray, “Failure time in heterogeneous systems,” Phys. Rev. Res. 1, 033047 (2019).
  * Moreira  _et al._ (2012) A. A. Moreira, C. L. N. Oliveira, A. Hansen, N. A. M. Araújo, H. J. Herrmann, and J. S. Andrade, “Fracturing highly disordered materials,” Phys. Rev. Lett. 109, 255701 (2012).
  * Roy  _et al._ (2015) Chandreyee Roy, Sumanta Kundu, and S. S. Manna, “Fiber bundle model with highly disordered breaking thresholds,” Phys. Rev. E 91, 032103 (2015).
  * Hope  _et al._ (2015) Sigmund Hope, Sumanta Kundu, Chandreyee Roy, Subhrangshu Manna, and Alex Hansen, “Network topology of the desert rose,” Frontiers in Physics 3 (2015), 10.3389/fphy.2015.00072.
  * George Batrouni and Hansen (1988) Ghassan George Batrouni and Alex Hansen, “Fourier acceleration of iterative processes in disordered systems,” Journal of Statistical Physics 52, 747–773 (1988).
  * Shekhawat  _et al._ (2013) Ashivni Shekhawat, Stefano Zapperi, and James P. Sethna, “From damage percolation to crack nucleation through finite size criticality,” Phys. Rev. Lett. 110, 185505 (2013).
  * Hansen and Schmittbuhl (2003) Alex Hansen and Jean Schmittbuhl, “Origin of the universal roughness exponent of brittle fracture surfaces:stress-weighted percolation in the damage zone,” Phys. Rev. Lett. 90, 045504 (2003).
  * Das and Biswas (2023) Soumyaditya Das and Soumyajyoti Biswas, “Critial scaling through gini-index,” arXiv:2211.01281 (2023).
  * Biswas and Chakrabarti (2020) Soumyajyoti Biswas and Bikas K. Chakrabarti, “Flory-like statistics of fracture in the fiber bundle model as obtained via kolmogorov dispersion for turbulence: A conjecture,” Phys. Rev. E 102, 012113 (2020).
  * Ramstad  _et al._ (2004) Thomas Ramstad, Jan Ø. H. Bakke, Johannes Bjelland, Torunn Stranden, and Alex Hansen, “Correlation length exponent in the three-dimensional fuse network,” Phys. Rev. E 70, 036123 (2004).


# Circling Back to Recurrent Models of Language

Gábor Melis
<EMAIL_ADDRESS>
DeepMind, London, UK; UCL, London, UK

###### Abstract

Just because some purely recurrent models suffer from being hard to optimize
and inefficient on today’s hardware, they are not necessarily bad models of
language. We demonstrate this by the extent to which these models can still be
improved by a combination of a slightly better recurrent cell, architecture,
objective, as well as optimization. In the process, we establish a new state
of the art for language modelling on small datasets and on enwik8 with dynamic
evaluation.

##  1 Introduction

Reliable model comparison is crucial for continued innovation in language
modelling. With so much attention on Transformers (Vaswani et al., 2017), the
development of purely recurrent models has ebbed away. It might be tempting to
claim that purely recurrent models are fundamentally worse models of language
than attention-based ones, but they seem to have an edge on small datasets,
while on larger datasets their lack of scalability on current hardware
(Hooker, 2020) also plays an important role in their evaluation. Due to being
unfashionable and slow, their fitness might be underestimated especially on
all but the smallest datasets with the exception of S4 (Gu et al., 2021), a
highly parallelizable, long-range model. The purpose of this work is simply to
apply more resources to advancing and evaluating recurrent models – despite,
and to compensate for, their computational ineffeciency – to better represent
their performance in future model comparisons.

Starting from the Mogrifier LSTM (Melis et al., 2019), a state-of-the-art
recurrent language model, we introduce several small changes, which involve
not only the design of the recurrent cell but the overall architecture, the
training objective, and optimization. On all datasets, they combine to a
significant, and in some cases large, effect.

##  2 Rewired LSTM

We propose a novel recurrent cell, a slightly tweaked version of the LSTM
(Hochreiter and Schmidhuber, 1997). In an LSTM, the updated state $\bm{c}$ and
the output $\bm{h}$ are computed as

$\displaystyle\bm{i}$
$\displaystyle=\sigma(\mathbf{W}^{ix}\bm{x}+\mathbf{W}^{ih}\bm{h}_{\textit{prev}}+\bm{b}^{i})$
$\displaystyle\bm{j}$
$\displaystyle=\tanh(\mathbf{W}^{jx}\bm{x}+\mathbf{W}^{jh}\bm{h}_{\textit{prev}}+\bm{b}^{j})$
$\displaystyle\bm{f}$
$\displaystyle=\sigma(\mathbf{W}^{fx}\bm{x}+\mathbf{W}^{fh}\bm{h}_{\textit{prev}}+\bm{b}^{f})$
$\displaystyle\bm{c}$
$\displaystyle=\bm{f}\odot\bm{c}_{\textit{prev}}+\bm{i}\odot\bm{j}$
$\displaystyle\bm{o}$
$\displaystyle=\sigma(\mathbf{W}^{ox}\bm{x}+\mathbf{W}^{oh}\bm{h}_{\textit{prev}}+\bm{b}^{o})$
$\displaystyle\bm{h}$ $\displaystyle=\bm{o}\odot\tanh(\bm{c}),$

where $\sigma$ is the logistic sigmoid function, $\odot$ is the elementwise
product, $\mathbf{W}^{**}$ and $\bm{b}^{*}$ are weight matrices and biases.
With $m$ and $n$ being the sizes of the input and the cell state, let the
$\operatorname*{LSTM}\colon\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}\times\mathbb{R}^{n}$
function refer to the cell computation given above:
$\operatorname*{LSTM}(\bm{c}_{\textit{prev}},\bm{h}_{\textit{prev}},\bm{x})=(\bm{c},\bm{h})$.
To accommodate the performance characteristics of parallel hardware, the
activations of $\bm{i}$, $\bm{j}$, $\bm{f}$, $\bm{o}$ are in practice often
computed by tiling the four $\mathbf{W}^{**}$ into one large matrix and
multiplying it with the concatenated input and recurrent state
$[\bm{x},\bm{h}_{\textit{prev}}]$.

Intuitively, in $\bm{c}=\bm{f}\odot\bm{c}_{\textit{prev}}+\bm{i}\odot\bm{j}$
the forget gate $\bm{f}$ and the “proposed update” $\bm{i}\odot\bm{j}$ depend
on each other, and computing them in parallel from the same inputs may be
partially redundant. On the flipside, the parametrization of the proposed
update as a product is more expressive than that of the forget gate, so the
overall cell update may lose some of the extra expressivity. Hence, it makes
sense to compute the forget gate from the proposed update instead. Note that
this reduces the opportunities for parallelization and makes cell updates
slower.

Further disregarding efficiency on today’s hardware, to reduce parameter count
and to encourage storing more information in the cell state $\bm{c}$, we
compute $\bm{o}$ from $\bm{c}$, dropping a potential bypass connection from
$\bm{x}$ to $\bm{h}$. Finally, to make exploding gradients less likely, we cap
$\bm{i}$ at $\bm{1-f}$ to ensure $|\bm{c}_{u}|\leqslant 1$ for all memory
units $u$. The update of our Rewired LSTM (RLSTM) cell takes the form

$\displaystyle\bm{i}$
$\displaystyle=\sigma(\mathbf{W}^{ix}\bm{x}+\mathbf{W}^{ih}\bm{h}_{\textit{prev}}+\bm{b}^{i})$
$\displaystyle\bm{j}$
$\displaystyle=\tanh(\mathbf{W}^{jx}\bm{x}+\mathbf{W}^{jh}\bm{h}_{\textit{prev}}+\bm{b}^{j})$
$\displaystyle\bm{f}$
$\displaystyle=\sigma(\mathbf{W}^{fu}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}i\odot
j}+\mathbf{W}^{fh}\bm{h}_{\textit{prev}}+\bm{b}^{f})$ $\displaystyle\bm{c}$
$\displaystyle=\bm{f}\odot\bm{c}_{\textit{prev}}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\min(\bm{i},1-\bm{f})}\odot\bm{j}$
$\displaystyle\bm{o}$
$\displaystyle=\sigma({\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathbf{W}^{oc}\bm{c}}+\bm{b}^{o})$
$\displaystyle\bm{h}$ $\displaystyle=\bm{o}\odot\tanh(\bm{c}),$

where the changes from the LSTM are highlighted in blue. Based on this altered
computation, we define the $\operatorname*{RLSTM}$ function similarly to the
$\operatorname*{LSTM}$ above.

##  3 Architecture

The way recurrent cells are combined can also be improved. Here, we simply opt
to use residual connections, where previous works have used stacked LSTMs
(Merity et al., 2017) or skip connections (Melis et al., 2019), that feed
directly into the final output. Along with this change, we apply dropout
(Hinton et al., 2012) to the cell output before it is added to the residual
branch.

To describe the overall architecture more formally, let us denote the
mogrification operation (Melis et al., 2019) with
$\operatorname*{mogrify}\colon\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}\times\mathbb{R}^{m}$:

$\displaystyle\operatorname*{mogrify}(\bm{h},\bm{x})$
$\displaystyle=\bm{h}^{2\lfloor r/2\rfloor},\bm{x}^{2\lfloor(r+1)/2\rfloor-1}$
$\displaystyle\bm{x}^{-1},\bm{h}^{0}$ $\displaystyle=\bm{x},\bm{h}$
$\displaystyle\bm{x}^{i}$
$\displaystyle=2\sigma(\mathbf{Q}^{i}\bm{h}^{i-1})\odot\bm{x}^{i-2},$
$\displaystyle\text{for odd i}\in[1\dots r]$ $\displaystyle\bm{h}^{i}$
$\displaystyle=2\sigma(\mathbf{R}^{i}\bm{x}^{i-1})\odot\bm{h}^{i-2},$
$\displaystyle\text{for even i}\in[1\dots r]$

where $r$ is a hyperparameter, and $\mathbf{Q}^{i}\in\mathbb{R}^{m\times
n},\mathbf{R}^{i}\in\mathbb{R}^{n\times m}$ are (possibly low-rank factorized)
weight matrices.

With time steps $t\in[1,2,\dots]$ and layers $l\in[1,L]$, from the vector of
token indices $\bm{w}$ in a fixed vocabulary, the probability distribution of
the next token $p(.|\bm{w}_{<t},\mathbf{M}_{t})$ is computed as follows.

$\displaystyle\bm{c}^{l}_{0},\bm{h}^{l}_{0}$ $\displaystyle=\bm{0},\bm{0}$
$\displaystyle\hat{\bm{x}}^{0}_{t}$
$\displaystyle=\operatorname*{onehot}(\bm{w}_{t})\mathbf{E}^{\textrm{in}}\odot\mathbf{M}^{\textrm{in}}_{t}$
$\displaystyle\bm{c}^{1}_{t},\bm{h}^{1}_{t}$
$\displaystyle=\textrm{[R]LSTM}(\bm{c}^{0}_{t-1},\operatorname*{mogrify}(\hat{\bm{h}}^{1}_{t-1},\hat{\bm{x}}^{0}_{t}))$
$\displaystyle\bm{c}^{l}_{t},\bm{h}^{l}_{t}$
$\displaystyle=\textrm{[R]LSTM}(\bm{c}^{l}_{t-1},\operatorname*{mogrify}(\hat{\bm{h}}^{l}_{t-1},\sum_{i=1}^{l-1}\hat{\bm{x}}^{i}_{t}))$
$\displaystyle(l>1)$ $\displaystyle\hat{\bm{x}}^{l}_{t}$
$\displaystyle=\bm{h}^{l}_{t}\odot\mathbf{M}^{\textrm{cell},l}_{t}$
$\displaystyle\hat{\bm{h}}^{l}_{t}$
$\displaystyle=\bm{h}^{l}_{t}\odot\mathbf{M}^{\textrm{state},l}$
$\displaystyle p(.|\bm{w}_{<t},\mathbf{M}_{t})$
$\displaystyle=\operatorname*{softmax}\Bigg{(}\bigg{(}\sum_{l=1}^{L}\hat{\bm{x}}^{l}_{t}\bigg{)}\odot\mathbf{M}^{\textrm{out}}_{t}\mathbf{E}^{\textrm{out}}+\bm{b}^{\textrm{out}}\Bigg{)},$

where we assume that $m=n$ to allow for a residual architecture without
projections. $\mathbf{E}^{\textrm{in}}$ and $\mathbf{E}^{\textrm{out}}$ are
the input and output embedding matrices. For word-based language modelling, we
set $\mathbf{E}^{\textrm{out}}$ to the transpose of $\mathbf{E}^{\textrm{in}}$
(Zoph and Le, 2016; Press and Wolf, 2016). $\mathbf{M}_{t}$ is the set of
individual input, cell, state and output dropout mask matrices
$\mathbf{M}^{\textrm{in}}_{t}$, $\mathbf{M}^{\textrm{cell},l}_{t}$,
$\mathbf{M}^{\textrm{state},l}$ and $\mathbf{M}^{\textrm{out}}_{t}$. Note that
for state dropout, we use the variational dropout of Gal and Ghahramani
(2016), thus $\mathbf{M}^{\textrm{state},l}$ does not depend on $t$. In
addition, when an RLSTM is used instead of an LSTM cell, we also apply the
state dropout mask to $\bm{c}$ in the calculation of
$\bm{o}=\sigma(\mathbf{W}^{oc}(\bm{c}\odot\mathbf{M}^{\textrm{state},l})+\bm{b}^{o})$
of the corresponding layers.

##  4 Objective

In the objective, we average model predictions over multiple dropout samples:

$\displaystyle\ln
p(\bm{w}_{t}|\bm{w}_{<t})=\ln\bigg{(}\frac{1}{D}\sum_{d=1}^{D}p(\bm{w}_{t}|\bm{w}_{<t},\mathbf{M}^{d}_{t})\bigg{)},$

where $D$ is the number of samples taken. As pointed out by Noh et al. (2017),
when dropout is interpreted as optimizing a variational lower bound on the log
likelihood (Gal and Ghahramani, 2016), this procedure is an instantiation of
Importance Weighted Autoencoders (Burda et al., 2015), which provide a bound
tighter than the single-sample ELBO.

##  5 Optimization

To increase the stability of optimization and allow slightly higher learning
rates, we use Rectified Adam (Liu et al., 2019). If training diverges, we
reset the weights and the optimization state to the previous best checkpoint
and multiply the learning rate by 0.9.

Merity et al. (2017) switch to averaging weights at a late stage of
optimization when the validation loss has not decreased for a while. A similar
procedure, called Stochastic Weight Averaging (SWA), was also proposed
(Izmailov et al., 2018), and corresponding theory was developed in Jain et al.
(2018) under the name of Tail Averaging. Here, we employ Two-Tailed Averaging
(Melis, 2022), which has no hyperparameters and provides a good approximation
to the optimal weight average at every step of optimization. Two-Tailed
Averaging (TTA) thus requires no tuning and is a much better fit with early
stopping, which is what we do on Enwik8 and Text8 (Hutter, 2012). While easier
to work with, TTA does not improve the final results over well-tuned Tail
Averaging, which is also used by our baseline, the Mogrifier LSTM.

LSTM and RLSTM forget gates are initialized with Chrono init (Tallec and
Ollivier, 2018) as $\bm{b}^{f}\sim\ln(\mathcal{U}(1,T_{\text{max}}-1))$, where
$T_{\text{max}}$ is a tuned hyperparameter.

Finally, we just train models longer where it is beneficial.

##  6 Dynamic Evaluation

Hinton and Plaut (1987) proposed fast weights, wherein parameters have a slow-
and a fast-changing component. Much later, Ba et al. (2016) showed a form of
attention to be an instantiation of fast weights. Similarly, some forms of
meta learning, e.g. few-shot adaptation with gradient updates, can be
interpreted as fast weights. To mimic this setting and gain a particularly
general fast weights implementation, it would be desirable to allow the model
weights to depend on the context as in $p(x_{i}|\theta(x_{<i}),x_{<i})$, where
the fast weights $\theta(x_{<i})$ are computed with a gradient-based updates
to the slow weights $\theta_{0}$. However, due to the practical difficulties
involved in training batches with per-example inner gradient updates, we
eschew fast weights at training time but not when performing evaluation, thus
we end up with what’s called dynamic evaluation (Krause et al., 2017). In this
work, we emphasize dynamic evaluation results as a proxy for the gradient-
based fast weights model.

As argued heuristically above, attention may be interpreted as a particular
way of adapting the weights to the context (outer product memory), like the
Mogrifier, which has an even more restricted form of update (scaling columns
of weight matrices). Thus, it is not surprising that Melis et al. (2019) found
that Transformers are better at adapting without changing their weights,
leaving less in-context signal for dynamic evaluation to pick up.

##  7 Experiments

###  7.1 Setup

We follow the experimental setup of our baseline (Melis et al., 2019). In the
following, we list only the most pertinent choices in our experimental setup;
everything else is the same as in the baseline. Note that the baseline’s and
our LSTM implementation already includes the capped input gate
($\min(1-\bm{f},\bm{i})$) of the RLSTM.

For the black-box hyperparameter tuner (Golovin et al., 2017), due to the
switch to a residual architecture, the baseline’s _inter_layer_dropout_
hyperparameter is replaced with _cell_output_dropout_ (see
$\mathbf{M}^{\textrm{cell},l}$ in Section 3). In addition, the top of the
Chrono init range $T_{\textit{max}}$ is a new hyperparameter from the
$[e^{2},e^{5}]$ range.

For word-level language modelling on Penn Treebank (Marcus et al., 1993) with
preprocessing by (Mikolov et al., 2010), we trained 2-layer models for about
400 epochs with 8 dropout samples. Batch size was 128, and we trained with a
BPTT (Werbos, 1990) window size of 70. Experiments on Wikitext-2 (Merity et
al., 2016) were conducted similarly, with the exception of training for 250
epochs and using 4 dropout samples.

For the character-based modelling task on Penn Treebank, we trained 2-layer
models for 100 epochs with 4 dropout samples, batch size 128, and BPTT window
size 200.

On Enwik8 and Text8 (Hutter, 2012), we trained 6-layer models for 200 epochs,
whereas the baseline model had 4 layers and was trained for only 29 epochs.
Batch size was 128, and the BPTT window size was set to 256. Increasing the
window size had no discernible effect in agreement with the findings of
Khandelwal et al. (2018). Due to the long time required to train these models,
the best hyperparameters were selected from a random pool of 60 candidates.
Since preliminary experiments indicated that the benefit of using multiple
dropout samples was less than 0.01 bpc, we refrained from using more than one
dropout sample to save time.

Model evaluation was performed with the standard, deterministic dropout; we
refrained from using the more expensive Monte-Carlo averaging (Gal and
Ghahramani, 2016). The optimal softmax temperature was selected at evaluation
time to maximize the validation log-likelihood (Melis et al., 2018). Finally,
we report results with and without dynamic evaluation (Krause et al., 2017).

###  7.2 Results

Due to the aforementioned loss of parallelization opportunities, we found that
the RLSTM was about 10%-30% slower than the LSTM. Still, the RLSTM retained a
significant advantage over the LSTM even when trained for the same wall clock
time.

Table 1: Word-level perplexities of near state-of-the-art models. Names for models with our new results are in bold. On Wikitext-2, the baseline employed Mixture of Softmaxes (Yang et al., 2017), but we found no benefit to that when used in conjunction with multiple dropout samples. |  |  | No Dyneval | Dyneval
---|---|---|---|---
|  |  | Val. | Test | Val. | Test
PTB EN |  Transformer-XL (Dai et al., 2019) | 24M | 56.7 | 54.5 |  |
Mogrifier LSTM (Melis et al., 2019) | 24M | 52.1 | 51.0 | 45.1 | 45.0
Mogrifier LSTM | 24M | 49.9 | 48.5 | 43.5 | 43.3
| Mogrifier RLSTM | 24M | 48.9 | 47.9 | 42.9 | 42.9
WT2 EN | Mogrifier LSTM MoS2 | 35M | 58.7 | 56.6 | 40.6 | 39.0
Mogrifier LSTM | 35M | 57.4 | 55.8 | 40.0 | 38.6
Mogrifier RLSTM | 35M | 56.7 | 55.0 | 39.3 | 38.0

On word-level language modelling (see Table 1), only training for half the
epochs, our Mogrifier LSTM outperformed the same in Melis et al. (2019). This
was mostly due to the quicker convergence when using multiple dropout samples
and, to a much smaller degree, to the initialization of the forget gate. Note
that on 2-layer models, that we used on this task, residual connections are
identical to skip connections, which are employed in Melis et al. (2019). On
top of these improvements, the RLSTM outperformed the LSTM by a small margin,
and we established a new state of the art on both datasets with and without
dynamic evaluation.

Table 2: Bits per character on character-based datasets of near state-of-the-art models. Names for models with our new results are in bold. The best test results with and without dynamic evaluation for a given dataset and model size are in bold unless there is a smaller model with a better result. |  |  | No Dyneval | Dyneval
---|---|---|---|---
|  |  | Val. | Test | Val. | Test
PTB EN |  Mogrifier LSTM (Melis et al., 2019) | 24M | 1.149 | 1.131 | 1.098 | 1.088
Mogrifier LSTM | 24M | 1.128 | 1.109 | 1.085 | 1.074
Mogrifier RLSTM | 24M | 1.115 | 1.096 | 1.073 | 1.061
Enwik8 EN |  Transformer-XL (d24) (Dai et al., 2019) | 277M |  | 0.993 |  | 0.940
Longformer (Beltagy et al., 2020) | 277M |  | 0.97 |  |
Transformer-XL (d18) (Dai et al., 2019) | 88M |  | 1.03 |  |
|  Longformer (Beltagy et al., 2020) | 102M |  | 0.99 |  |
|  Mogrifier LSTM (Melis et al., 2019) | 96M | 1.110 | 1.122 | 1.009 | 0.988
| Mogrifier LSTM | 96M | 1.057 | 1.073 | 0.963 | 0.946
| Mogrifier RLSTM | 96M | 1.028 | 1.042 | 0.952 | 0.935
|  Transformer-XL (d12) (Dai et al., 2019) | 41M |  | 1.06 |  | 1.01
|  Longformer (Beltagy et al., 2020) | 41M | 1.02 | 1.00 |  |
|  Mogrifier LSTM (Melis et al., 2019) | 48M | 1.135 | 1.146 | 1.035 | 1.012
| Mogrifier LSTM | 48M | 1.083 | 1.094 | 0.988 | 0.970
| Mogrifier RLSTM | 48M | 1.060 | 1.071 | 0.986 | 0.968
Text8 EN |  Transformer-XL (d24) (Dai et al., 2019) | 277M |  | 1.08 |  | 1.038
Mogrifier LSTM | 96M | 1.033 | 1.106 | 0.977 | 1.047
Mogrifier RLSTM | 96M | 1.022 | 1.096 | 0.975 | 1.044
|  Longformer (Beltagy et al., 2020) | 41M | 1.04 | 1.10 |  |
| Mogrifier LSTM | 48M | 1.063 | 1.140 | 1.005 | 1.075
| Mogrifier RLSTM | 48M | 1.044 | 1.119 | 0.998 | 1.068

Table 2 shows our results on character-based language modelling. On Penn
Treebank, again only training for half the epochs, we significantly boosted
the results of the state-of-the-art baseline model by using multiple dropout
samples and to a smaller degree by tuning the forget gate initialization. Our
results were then further improved by switching to the RLSTM.

Our results on Enwik8 and Text8 were boosted greatly. However, some corners
were cut due the high computational cost, thus the results can likely be
improved further e.g. with proper tuning (instead of random sampling), by
training even longer, and by using multiple dropout samples. We heuristically
estimate a 0.015-0.03 bpc suboptimality due to these factors only.

Taking advantage of TTA’s online estimates and the fact that hyperparameters
were selected randomly, we can estimate the contribution of training longer
(200 vs 29 epochs). For example, in the 48M parameter setting on Enwik8, we
observed that training longer lowered the validation bpc by about 0.035,
leaving another 0.016 bpc for the contributions of residual connections and
the forget gate initialization. Finally, the RLSTM outperformed the LSTM by
0.022 bpc. Similar relative contributions were observed in all other settings
on Enwik8 and Text8.

Some of the changes we made increased primarily the stability of training and
the efficiency of hyperparameter tuning without discernible effect on the
final results. Using Rectified Adam and restarting from a previous checkpoint
on divergence greatly reduced the number of failed runs, and Two-Tailed
Averaging (TTA) made tuning easier by removing a hyperparameter while opening
the door for early stopping due to its online nature.

We found that in the early stages using $D$ dropout samples allowed
optimization to make more progress per step but less than the progress with a
single sample and $D$ times the number of steps. On datasets which are on the
smaller side and where dropout rates are high, in the late stages of
optimization, the multi-sample objective proved better in terms of validation
perplexity. However, that was not the case on Enwik8 and Text8, possibly due
to being bottlenecked by other issues such as the length of optimization.

##  8 Conclusion

We strengthened purely recurrent models’ language modelling results on a
varied collection of datasets using various techniques and a slightly novel
recurrent cell. These new baselines better represent the models’ ability but
do not necessarily allow for a fair comparison to previous results. In
particular, while we cited and listed results of the best transformer-based
models, we did not evaluate them ourselves using the same methodology. This is
especially obvious where previous works did not report dynamic evaluation
results. The comparisons we made to recurrent models were also rather
targeted: we focussed solely on the state-of-the-art purely recurrent model,
the Mogrifier LSTM, and ignored both the less performant models and those
combined with attention (Merity, 2019; Lei, 2021). Nevertheless, the
combination of our strong results and little novelty highlight the extent to
which minor details and computational efficiency on current hardware affect
model comparisons.

With that in mind, the primary contribution of this work is to apply more
resources to designing and evaluating recurrent models and to provide stronger
baselines for model comparisons. Second, the set of improvements to models of
the recurrent cell, overall architecture, training objective, and optimization
that we employed or proposed might inform future practice and research.

## References

  * Ba et al. (2016) Jimmy Ba, Geoffrey E Hinton, Volodymyr Mnih, Joel Z Leibo, and Catalin Ionescu.  Using fast weights to attend to the recent past.  _Advances in neural information processing systems_ , 29, 2016.
  * Beltagy et al. (2020) Iz Beltagy, Matthew E Peters, and Arman Cohan.  Longformer: The long-document transformer.  _arXiv preprint arXiv:2004.05150_ , 2020.
  * Burda et al. (2015) Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov.  Importance weighted autoencoders.  _arXiv preprint arXiv:1509.00519_ , 2015.
  * Dai et al. (2019) Zihang Dai, Zhilin Yang, Yiming Yang, William W Cohen, Jaime Carbonell, Quoc V Le, and Ruslan Salakhutdinov.  Transformer-XL: Attentive language models beyond a fixed-length context.  _arXiv preprint arXiv:1901.02860_ , 2019.
  * Gal and Ghahramani (2016) Yarin Gal and Zoubin Ghahramani.  A theoretically grounded application of dropout in recurrent neural networks.  In _Advances in Neural Information Processing Systems_ , pages 1019–1027, 2016.
  * Golovin et al. (2017) Daniel Golovin, Benjamin Solnik, Subhodeep Moitra, Greg Kochanski, John Karro, and D Sculley.  Google Vizier: A service for black-box optimization.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1487–1495. ACM, 2017.
  * Gu et al. (2021) Albert Gu, Karan Goel, and Christopher Ré.  Efficiently modeling long sequences with structured state spaces.  _arXiv preprint arXiv:2111.00396_ , 2021.
  * Hinton and Plaut (1987) Geoffrey E Hinton and David C Plaut.  Using fast weights to deblur old memories.  In _Proceedings of the 9th annual conference of the cognitive science society_ , pages 177–186, 1987.
  * Hinton et al. (2012) Geoffrey E Hinton, Nitish Srivastava, Alex Krizhevsky, Ilya Sutskever, and Ruslan R Salakhutdinov.  Improving neural networks by preventing co-adaptation of feature detectors.  _arXiv preprint arXiv:1207.0580_ , 2012.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber.  LSTM can solve hard long time lag problems.  In _Advances in neural information processing systems_ , pages 473–479, 1997.
  * Hooker (2020) Sara Hooker.  The hardware lottery.  _CoRR_ , abs/2009.06489, 2020.  URL https://arxiv.org/abs/2009.06489.
  * Hutter (2012) Marcus Hutter.  The human knowledge compression contest.  _URL http://prize. hutter1. net_ , 6, 2012.
  * Izmailov et al. (2018) Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov, and Andrew Gordon Wilson.  Averaging weights leads to wider optima and better generalization.  _arXiv preprint arXiv:1803.05407_ , 2018.
  * Jain et al. (2018) Prateek Jain, Sham Kakade, Rahul Kidambi, Praneeth Netrapalli, and Aaron Sidford.  Parallelizing stochastic gradient descent for least squares regression: mini-batching, averaging, and model misspecification.  _Journal of Machine Learning Research_ , 18, 2018.
  * Khandelwal et al. (2018) Urvashi Khandelwal, He He, Peng Qi, and Dan Jurafsky.  Sharp nearby, fuzzy far away: How neural language models use context.  _arXiv preprint arXiv:1805.04623_ , 2018.
  * Krause et al. (2017) Ben Krause, Emmanuel Kahembwe, Iain Murray, and Steve Renals.  Dynamic evaluation of neural sequence models.  _arXiv preprint arXiv:1709.07432_ , 2017.
  * Lei (2021) Tao Lei.  When attention meets fast recurrence: Training language models with reduced compute.  _CoRR_ , abs/2102.12459, 2021.  URL https://arxiv.org/abs/2102.12459.
  * Liu et al. (2019) Liyuan Liu, Haoming Jiang, Pengcheng He, Weizhu Chen, Xiaodong Liu, Jianfeng Gao, and Jiawei Han.  On the variance of the adaptive learning rate and beyond.  _arXiv preprint arXiv:1908.03265_ , 2019.
  * Marcus et al. (1993) Mitchell P Marcus, Mary Ann Marcinkiewicz, and Beatrice Santorini.  Building a large annotated corpus of english: The Penn treebank.  _Computational linguistics_ , 19(2):313–330, 1993\.
  * Melis et al. (2018) Gábor Melis, Charles Blundell, Tomáš Kočiskỳ, Karl Moritz Hermann, Chris Dyer, and Phil Blunsom.  Pushing the bounds of dropout.  _arXiv preprint arXiv:1805.09208_ , 2018.
  * Melis et al. (2019) Gábor Melis, Tomáš Kočiskỳ, and Phil Blunsom.  Mogrifier LSTM.  _arXiv preprint arXiv:1909.01792_ , 2019.
  * Melis (2022) Gábor Melis.  Two-tailed averaging: Anytime adaptive once-in-a-while optimal iterate averaging for stochastic optimization, 2022.  URL https://arxiv.org/abs/2209.12581.
  * Merity (2019) Stephen Merity.  Single headed attention rnn: Stop thinking with your head.  _arXiv preprint arXiv:1911.11423_ , 2019.
  * Merity et al. (2016) Stephen Merity, Caiming Xiong, James Bradbury, and Richard Socher.  Pointer sentinel mixture models.  _CoRR_ , abs/1609.07843, 2016.  URL http://arxiv.org/abs/1609.07843.
  * Merity et al. (2017) Stephen Merity, Nitish Shirish Keskar, and Richard Socher.  Regularizing and optimizing LSTM language models.  _arXiv preprint arXiv:1708.02182_ , 2017.
  * Mikolov et al. (2010) Tomas Mikolov, Martin Karafiát, Lukas Burget, Jan Cernockỳ, and Sanjeev Khudanpur.  Recurrent neural network based language model.  In _Interspeech_ , volume 2, page 3, 2010.
  * Noh et al. (2017) Hyeonwoo Noh, Tackgeun You, Jonghwan Mun, and Bohyung Han.  Regularizing deep neural networks by noise: Its interpretation and optimization.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Press and Wolf (2016) Ofir Press and Lior Wolf.  Using the output embedding to improve language models.  _CoRR_ , abs/1608.05859, 2016.  URL http://arxiv.org/abs/1608.05859.
  * Tallec and Ollivier (2018) Corentin Tallec and Yann Ollivier.  Can recurrent neural networks warp time?  _arXiv preprint arXiv:1804.11188_ , 2018.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * 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.
  * Yang et al. (2017) Zichao Yang, Zhiting Hu, Ruslan Salakhutdinov, and Taylor Berg-Kirkpatrick.  Improved variational autoencoders for text modeling using dilated convolutions.  In _International conference on machine learning_ , pages 3881–3890. PMLR, 2017.
  * Zoph and Le (2016) Barret Zoph and Quoc V. Le.  Neural architecture search with reinforcement learning.  _CoRR_ , abs/1611.01578, 2016.  URL http://arxiv.org/abs/1611.01578.