# Error-Correcting Neural Sequence Prediction

James O’ Neill and Danushka Bollegala
Department of Computer Science, University of Liverpool
Liverpool, L69 3BX
England
{james.o-neill<EMAIL_ADDRESS>

###### Abstract

We propose a novel neural sequence prediction method based on error-correcting
output codes that avoids exact softmax normalization and allows for a tradeoff
between speed and performance. Instead of minimizing measures between the
predicted probability distribution and true distribution, we use error-
correcting codes to represent both predictions and outputs. Secondly, we
propose multiple ways to improve accuracy and convergence rates by maximizing
the separability between codes that correspond to classes proportional to word
embedding similarities.

Lastly, we introduce our main contribution called Latent Variable Mixture
Sampling, a technique that is used to mitigate exposure bias, which can be
integrated into training latent variable-based neural sequence predictors such
as ECOC. This involves mixing the latent codes of past predictions and past
targets in one of two ways: (1) according to a predefined sampling schedule or
(2) a differentiable sampling procedure whereby the mixing probability is
learned throughout training by replacing the greedy argmax operation with a
smooth approximation. ECOC-NSP leads to consistent improvements on language
modelling datasets and the proposed Latent Variable mixture sampling methods
are found to perform well for text generation tasks such as image captioning.

##  1 Introduction

Sequence modelling (SM) is a fundamental task in natural language which
requires a parametric model to generate tokens given past tokens. SM underlies
many types of structured prediction tasks in natural language, such as
Language Modelling (LM), Sequence Tagging (e.g Named Entity Recognition,
Constituency/Dependency Parsing) and Text Generation (e.g Image Captioning and
Machine Answering (Sutskever et al., 2014)). The goal is to learn a joint
probability distribution for a sequence of length $T$ containing words from a
vocabulary $\mathcal{V}$. This distribution can be decomposed into the
conditional distributions of current tokens given past tokens using the chain
rule shown in Equation 1. In Neural Sequence Prediction (NSP), a Recurrent
Neural Network (RNN) $f_{\theta}(\cdot)$ parameterized by $\theta$ is used to
encode the information at each timestep $t$ into a hidden state vector
$h^{l}_{t}$ which is followed by a decoder $z^{l}_{t}=h^{l}_{t}W^{l}+b^{l}$
and a normalization function $\phi(z_{t}^{l})$ which forms a probability
distribution $\hat{p}_{\theta}(y_{t}|x_{t},h_{t-1})$, $\;\forall
t\in[0,1,..T]$.

$P(w_{1},...,w_{T})=\prod^{T}_{t=1}P(w_{t}|w_{t−1},...,w_{1})$ (1)

However, training can be slow when $|\mathcal{V}|$ is large while also leaving
a large memory footprint for the respective input embedding matrices.
Conversely, in cases where the decoder is limited by an information
bottleneck, the opposite is required where more degrees of freedom are
necessary to alleviate information loss in the decoder bottleneck. Both
scenarios correspond to a trade-off between computation complexity and out-of-
sample performance. Hence, we require that a newly proposed model has the
property that the decoder can be easily configured to deal with this trade-off
in a principled way. Lastly, standard supervised learning (self-supervised for
sequence prediction) assumes inputs are i.i.d. However, in sequence
prediction, the model has to rely on its own predictions at test time, instead
of past targets that are used as input at training time. This difference is
known as exposure bias and can lead to errors compounding along a generated
sequence. This approach to sequence prediction is also known as teacher
forcing where the teacher provides all targets at training time. We also
require that exposure bias is addressed while dealing with the aforementioned
challenges related to computation and performance trade-offs in the decoder.
Hence, we propose an error-correcting output code (ECOC) based NSP (ECOC-NSP)
that address this desideratum. We show that when given sufficient error codes
($|\mathcal{V}|\gg|c|\gg\log_{2}(|\mathcal{V}|)$) while the codeword
dimensionality $|c|<|\mathcal{V}|$, accuracy can be maintained compared to
using the full softmax. Lastly, we show that this latent variable-based NSP
approach can be extended to mitigate the aforementioned problem of compounding
errors by using Latent Variable Mixture Sampling (LVMS). LVMS in an ECOC-NSP
model also outperforms an equivalent Hierarchical Softmax-based NSP that uses
Scheduled Sampling (Bengio et al., 2015) and other closely related baselines.
This is the first report of mitigating compounding errors when approximating
the posterior in recurrent neural networks (RNNs).

##### Contributions

Our main contributions are summarized as the following:

  1. 1.

An error-correcting output coded neural language model that requires less
parameters than its softmax-based sequence modelling counterpart given
sufficient separability between classes via error-checks.

  2. 2.

An embedding cosine similarity rank ordered codebook that leads to well-
separated codewords, where the number of error-correcting codes assigned to a
token is proportional to the cosine similarity between the tokens
corresponding pretrained word embedding and the most frequent tokens word
embedding.

  3. 3.

A Latent-Mixture Sampling method to mitigate exposure bias in Latent Variable
models. This is then extended to Differentiable Latent Variable Mixture
Sampling that uses the Gumbel-Softmax so that discrete categorical variables
can be backpropogated through. This performs comparably to other sampling-
based approaches.

  4. 4.

Novel baselines such as Scheduled Sampling Hierarchical Softmax (SS-HS) and
Scheduled Sampling Adaptive Softmax (SS-AS), are introduced in the evaluation
of our proposed ECOC method. This applies SS to closely related softmax
approximations.

##  2 Background

###  2.1 Error-Correcting Codes

Error-Correcting Codes (Hamming, 1950) originate from seminal work in solid-
state electronics around the time of the first digital computer. Later, binary
codes were introduced in the context of artificial intelligence via the
NETtalk system (Sejnowski and Rosenberg, 1987), where each class index is
represented by a binary code $C$ and a predicted code as $\hat{C}$ from some
parametric model $f_{\theta}(\cdot)$. When $|\mathcal{V}|\neq 2^{n}$, the
remaining codes are used as error-correction bits $k=|\mathcal{V}|-2^{n}$.
This tolerance can be used to account for the information loss due to the
sample size by measuring the distance (e.g Hamming) between the predicted
codeword and the true codeword with $d$ error-correction bits. If the minimum
distance between codewords is $d$ then at least $(d-1)/2$ bits can be
corrected for and hence, if the Hamming distance $d\leq(d-1)/2$ we will still
retrieve the correct codeword. In contrast to using one bit per $k$ classes in
standard multi-class classification, error-correction cannot be achieved. Both
error-correction and class bits make up the codebook $\mathcal{C}$.

##### Why Latent Codes for Neural Sequence Prediction?

Targets are represented as 1-hot vectors (i.e. kronecker delta) in standard
training of neural sequence predictors, treated as a 1-vs-rest multi-class
classification. This approach can be considered a special case of ECOC
classification where the codebook $\mathcal{C}$ with $n$ classes is
represented by an identity $\mathbb{I}^{n\times n}$. ECOC classification is
well suited over explicitly using observed variables when the output space is
structured. Flat-classification (1-vs-rest) ignores the dependencies in the
outputs, in comparison to using latent codes that share some common latent
variables between associative words. For example, in the former case, if we
train a model that only observes the word silver in a sequence …silver car..
and then at test-time observes silver-back, because there is high association
between silver and car, the model is more likely to predict car instead of
gorilla. ECOC is less prone to such mistakes because although a/some bit/s may
be different between the latent codes for car and gorilla, the potential
misclassifications can be re-corrected with the error-correcting bits. In
other words, latent coding can reduce the variance of each individual
classifier and has some tolerance to mistakes induced by sparse transitions,
proportional to the number of error-checks used.

###  2.2 Methods for Softmax Approximation

Goodman (2001); Morin and Bengio (2005) proposed a Hierarchical Softmax (HS)
that outputs short codes representing marginals of the conditional
distribution, where the product of these marginals along a path in the tree
approximate the conditional distribution. This speeds up training proportional
to the traversed binary tree path lengths, where intermediate nodes assign
relative probabilities of child nodes. Defining a good tree structure improves
performance since semantically similar words have a shorter path and therefore
similar representations are learned for similar words. HS and ECOC are similar
insofar as they both can be interpreted as approximating the posterior as a
product of marginal probabilities using short codes. However, ECOC tolerates
errors in code prediction. This becomes more important as the code length (i.e
the tree depth in HS) grows since the likelihood of mistakes becomes higher.
Hence, ECOC can be considered a flexible tradeoff between the strictness of HS
and the full softmax. The Differentiated Softmax (DS) uses a sparse linear
block of weights for the decoder where a set of partitions are made according
to the unigram distribution, where the number of weights are assigned
proportional to the term frequency. This is intuitive since rare words require
less degrees of freedom to account for the little amount of contexts in which
they appear, compared to common words. The optimal spacing of error-checking
bits between codewords corresponding to tokens is similar to - varying
branching sizes for different areas of the tree depending on the unigram
frequency or allocating weights proportional to the frequency akin to
Differentiable Softmax (Chen et al., 2015). We also consider achieving the
spacing via term frequency, but also use our proposed rank ordered word
embedding cosine similarities as mentioned in Equation 3.1.

The Adaptive Softmax (AS) Grave et al. (2016) provide an approximate
hierarchical model that directly accounts for the computation time of matrix
multiplications. AS results in 2x-10x speedups when compared to the standard
softmax, dependent on the size of the corpus and vocabulary. They find on
sufficiently large corpora (Text8, Europarl and 1-Billion datasets), accuracy
is maintained while reducing the computation time.

###  2.3 Recent Applications of Latent Codes

Shu and Nakayama (2017) recently used compositional codes for word embeddings
to cut down on memory requirements in mobile devices. Instead of using binary
coding, they achieve word embedding compression using multi-codebook
quantization. Each bit $c\in C$ comprises of a discrete code (0-9) and
therefore at minimum $\log_{10}(k)$ bits are required. They also propose to
use the Gumbel-Softmax trick but for the purposes of learning the discrete
codes. The performance was maintaned for sentiment analysis and machine
translation with 94% and 98% respective compression rates. Shi and Yu (2018)
propose a product quantizatioon structured embedding that reduces memory by
10-20 times the number of parameters, while maintaining performance. This
involves slicing the embedding tensor into groups which are then quantized
independently. Oda et al. (2017) also consider binary code for neural machine
translation. However, their approach still required binary codes to be used in
conjunction with the softmax and showed a performance degradation when only
using binary codes. Here, we show that, when given enough bits, the model is
competitive against the full softmax, and in some cases outperforming.
Moreover, we introduce novel ways to mitigate exposure bias in these models,
and Latent Variable based models alike.

##  3 Methodology

###  3.1 Codebook Construction

A challenging aspect of assigning codewords is ordering the codes so that even
if incorrect predictions are made, that the codeword is at least semantically
closer to that of the codewords that are less related, while ensuring good
separation between codes. Additionally, we have to consider the amount of
error-checking bits to use. In theory, $\log_{2}(k)/k$ is sufficient to
account for all $k$ classes. However, this alone can lead a degradation in
performance. Hence, we also consider a large amount of error-checking bits. In
this case, the error-checking bits can account for more mistakes given by
other classes, which may be correlated. In contrast, using probability
distributions naturally accounts for these correlations, as the mass needs to
shift relative to the activation of each output. This is particularly
important for language modelling and text generation because of the high-
dimensionality of the output. The most naive way to create the codebook is to
assign binary codes to each word in random order. However, it is preferable to
assign similar codes to $w\in\mathcal{V}$ that are semantically similar while
maximizing the Hamming distance between codes where leftover error codes
separate class codes.

####  3.1.1 Codebook Arrangement

A fundamental challenge in creating the codebook $\mathcal{C}$ is in how
error-codes are distributed between codes that maximize the separability
between codewords that are more likely to be interchangeably and incorrectly
predicted. This is related to the second challenge of choosing the
dimensionality of $C$. The latter is dependent on the size of the corpus, and
in some cases might only require $|\log_{2}(\mathcal{V})|\leq
d\leq|\mathcal{V}|$ bits to represent all classes with leftover error-checking
bits. These two decisions correspond to a tradeoff between computational
complexity and accuracy of our neural language model, akin to tree
expressitivity in the Hierarhcial Softmax to using the Full Softmax. Below we
describe a semantically motivated method to achieve well-separated codewords,
followed by a guide on how to choose codebook dimensionality
$d_{\mathcal{C}}$.

##### Embedding Simlarity-Based Codebooks

Previous work on ECOC has focused on theories as to why randomly generated
codes lead to good row and column separation (Berger, 1999). However, this
assumes that class labels are conditionally independent and therefore it does
not apply well for sequence modelling where the output space is loosely
structured. To address this, we propose to reorder $C\in\mathcal{C}$ such that
Hamming distance between any two codewords is proportional to the embedding
similarity. Moreover, separating codewords by semantic similarity can be
achieved by placing the amount of error-checking bits proportional to rank
ordered similarity for a chosen query word embedding. A codebook ordered by
pretrained word embedding similarities for $w_{*}$ is denoted as
$\mathcal{C}_{\Lambda_{w_{*}}}$. The similarity scores between embeddings is
given as $\mathcal{F}(\mathcal{M}_{*},M_{i})\;\forall i$ is used reorder
$M\to\mathcal{M}^{{}^{\prime}}$. In our experiments we use pretrained
GoogleNews skipgram embeddings111available here:
https://code.google.com/archive/p/word2vec/.

Good separation is achieved when codes are separated proportional to the
cosine similarity between embeddings of the most frequent word
$w_{*}\in\mathcal{V}$ and the remaining words $w^{\prime}$. Therefore, words
with high similarity have corresponding codes that are closer in Hamming
distance $\text{H}(\cdot,\cdot)$ in $\mathcal{C}$. This ensures that even when
codes are correlated, that incorrect latent predictions are at least more
likely to correspond to semantically related words. We are not guaranteed that
codes close in Hamming distance are closer in a semantic sense in the random
case. Therefore, we can instead consider computing $\mathcal{M}^{{}^{\prime}}$
using ordered similarities of word embedding similarities where the function
$\mathcal{F}(\cdot,\cdot)$ computes the cosine similarity for any two words.
If we apply the same combinatorial analysis to codebook separation when the
codes are assigned proportional to semantic similarity, we can define the good
separation as the following:

######  Definition 1.

Given $k$ redundant codewords $C_{k}$, we require an assignment that leads to
a strongly separated $\mathcal{C}$. Let $\delta(\cdot,\cdot)$ denote a
function that assigns $C_{i}^{k}$ error-checking codewords assigned to the
$i^{th}$ class codeword and
$\delta(C_{*},C_{i})\propto\mathcal{F}(\mathcal{M}_{*},\mathcal{M}_{i})\;\forall
k$.

In practice $\delta(\cdot,\cdot)$ normalizes the resulting embedding
similarities $S=\mathcal{F}(M_{*},M)$ using a normalization function cumsum
$\big{(}\phi(S)\times|\mathcal{C}|\big{)}$ to assign the intervals between
adjacent codeword spans. This assigns greater distance between words that are
more similar to $w*$, and less error-checking codewords to rare words that
have distant neighbouring words in the embedding space.

###  3.2 Latent Variable Mixture Sampling

To mitigate exposure bias for latent variable-based sequence modelling we
propose a sampling strategy that interpolates between predicted and target
codewords. We refer to this as Latent Variable Mixture Sampling (LVMS) and its
application to ECOC as Codeword Mixture Sampling (CMS).

##### Curriculum-Based Latent Variable Mixture Sampling

Figure 1: Curriculum Mixture Sampling

In Curriculum-Based Latent Variable Mixture Sampling (CLVMS), the mixture
probability is $p_{c}=0\;\forall c\in C$ at epoch $\epsilon=0$ and throughout
training the probability monotonically increases $p_{c}=\delta_{c}\;\forall
c\in C$ where $\delta_{c}$ is the threshold for the $c$ th bit after
$\epsilon$ epochs. A Bernoulli sample
$\tilde{C}=\mathbb{B}(\hat{C}_{c},C_{c})\;\forall c\in[0,C]$ is carried out
for $t\in T$ in each minibatch. The probabilities per dimension $p_{c}$ are
independent of keeping a prediction $\hat{y}_{t-1,c}$ instead of the $c$ th
bit in the target codeword $y_{(t-1,c)}$ at timestep $t$-$1$. The reason for
having individual mixture probabilities per bits is because when we consider a
default order in $\mathcal{C}$, this results in tokens being assigned
codewords ranked by frequency. Therefore, the leftmost bit predictions are
more significant than bit errors near the beginning (e.g $2^{0}=1$ only 1 bit
difference). We report results when using a sigmoidal schedule as shown in
Equation 2 where $\tau_{max}$ represents the temperature at the last epoch and
$\delta$ is a scaling factor controlling the slope.

$[\hat{y}_{t-1},y_{t-1}]\sim\frac{\tau_{max}}{1+\exp(-\epsilon/\delta)},\;\forall\epsilon\in[-\text{N}/2,\text{N}/2]$
(2)

This is different to scheduled sampling since we are not just choosing between
the ground truth and prediction, instead we are mixing factors of the
predicted factored distribution and the target factored distribution that
represents the posterior, as illustrated in Figure 1. Here, the color strength
illustrates the activation between $[0,1]$. Figure 1 also demonstrates how CMS
is used in a HS without the additional error checks, in which case each
activation corresponding to a node in the tree.

##### Latent Soft-Mixture Sampling

In standard CMS, we pass token index $w_{t}$ which is converted to an input
embedding $e_{w}$ based on the most probable bits predictions at the last time
step $\operatorname*{arg\,max}_{\theta}p(y_{t-1}|x_{t-1};\theta)$. We can
instead replace the argmax operator with a soft argmax that uses a weighted
average of embeddings $e\in E$ where weights are assigned from the previous
predicted output via the softmax normalization $\phi(x_{t-1},\tau)$, where
$\tau$ controls the kurtosis of the probability distribution ($\tau\to 0$
tends to argmax) in Equation 3.

$x_{t}=\sum_{w\in\mathcal{V}}e_{w}\Big{(}\frac{\exp(h^{T}_{w}\theta/\tau)}{\sum_{w\in\mathcal{V}}\exp(h^{T}_{w}\theta/\tau)}\Big{)}$
(3)

In the ECOC-NSP, we consider binary codewords and therefore choose the top $k$
least probable bits to flip according to the curriculum schedule. Hence, this
results in $k$ codewords where each $\hat{C}$ has at least hamming distance
$\text{H}(\hat{C},C)=1$ $(2^{0})$. Concretely, this is a soft interpolation
between past targets and a weighted sum of the $k$ most probable codewords
$\hat{C}_{K}=\operatorname*{arg\,max}_{k}\big{(}\sigma(h^{T}_{w}W)\big{)}$
such that $x_{t}=\mathbb{B}_{K}\Big{(}C,\sum_{k}^{K}\phi(\hat{C}_{k})\Big{)}$
where $B_{K}$ samples one or the other for each $k$th dimension of $C$.

###  3.3 Differentiable Latent Variable Sampling

The previous curriculum strategies disregard where the errors originate from.
Instead, they interpolate between model predictions of latent variables
$\hat{Y}$ and targets $Y$ in a way that does not distinguish between cascading
errors and localized errors. This means that it only recorrects errors after
they are made instead of directly correcting for the origin of the errors.
Maddison et al. (2016) showed that such operations can be approximated by
using a continuous relaxation using the reparameterization trick, also known
as the Concrete Distribution. By applying such relaxation it allows us to
sample from the distribution across codes while allowing for a fully
differentiable objective, similar to recent work (Goyal et al., 2017). We
extend this to mixture sampling by replacing the argmax operation with the
Concrete distribution to allow the gradients to be adjusted at points where
prior predictions changed value throughout training. This not only identifies
at which time-step the error occurs, but what latent variables (i.e. output
codes) had the most influence in generating the error. This is partially
motivated by the finding that in the latent variable formulation of simple
logistic regression models, the latent variable errors form a Gumbel
distribution. Hence, we sample latent codes inversely proportional to the
errors from a Gumbel distribution.

Gumbel-Softmax Similarly, instead of passing the most likely predicted word
$\hat{y}^{w*}_{t-1}$, we can instead sample from
$\hat{y}_{t-1}\sim\phi(h_{t-1},w)$ and then pass this index as $\hat{x}_{t}$.
This is an alternative to always acting greedily and allow the model to seek
other likely actions. However, to compute derivatives through samples from the
softmax, we need avoid discontinuities, such as the argmax operation. The
Gumbel-Softmax (Maddison et al., 2016; Jang et al., 2016) allows us to sample
and differentiate through the softmax by providing a continuous relaxation
that results in probabilities instead of a step function (i.e. argmax). As
shown in Equation 4, for each componentwise Gumbel noise $k\in[1..,n]$ for
latent variable given by $h^{T}\theta$, we find $k$ that maximizes
$\log\alpha_{k}-\log(-\log U_{k})$ and then set $D_{k}=1$ and $D{\neg k}=0$,
where $U_{k}\sim\text{Uniform}(0,1)$ and $\alpha_{k}$ is drawn from a discrete
distribution $D\sim\text{Discrete}(\alpha)$.

$\displaystyle\hat{p}(y_{t}|x_{t};\theta)=\frac{\exp((\log\alpha_{k}+G_{k})/\tau)}{\sum_{i=1}^{n}\exp((\log\alpha_{i}+G_{i})/\tau)}$
(4)

For ECOC, we instead consider Bernoulli random variables which for the
Concrete distribution can be expressed by means of two arbitrary Gumbel
distributions $G_{1}$ and $G_{2}$. The difference between $G_{1}$ and $G_{2}$
follows a Logistic distribution and so $G_{1}-G_{2}\sim\text{Logistic}$ and is
sampled as $G_{1}-G_{2}\equiv\log U-\log(1-U)$. Hence, if
$\alpha=\alpha_{1}/\alpha_{2}$, then
$P(D_{1}=1)=P(G_{1}+\log\alpha_{1}>G_{2}+\log\alpha_{2})=P(\log
U−\log(1-U)+\log\alpha>0)$. For a step function $\mathcal{H}$,
$D_{1}\equiv\mathcal{H}(\log\alpha+\log U-\log(1-U))$, corresponding to the
Gumbel Max-Trick (Jang et al., 2016).
Sampling a Binary Concrete random variable involves sampling $Z$, sample
$L\sim\text{Logistic}$ and set $Z$ as shown in Equation 5, where
$\alpha,\tau\in(0,\infty)$ and $Z\in(0,1)$. This Binary Concrete distribution
is henceforth denoted as $\text{BinConcrete}(\cdot,\cdot)$ with location
$\alpha$ and temperature $\tau$. In the forward pass the probability $Z$ is
used to compute the approximate posterior, unlike the one-hot vectors used in
straight-through estimation.

$Z\equiv\frac{1}{1+\exp\big{(}−(\log\alpha+L)/\tau\big{)}}$ (5)

This is used for ECOC and other latent variable-based models, such as
Hierarchical Sampling, to propogate through past decisions and make corrective
updates that backpropogate to where errors originated from along the sequence.
Hence, we also carry out experiments with BinConcrete (Equation 5) and Gumbel-
Softmax( Equation 4) for HS and ECOC respectively. The temperature $\tau$ can
be kept static, annealed according to a schedule or learned during training,
in the latter case this is equivalent to entropy regularization (Grandvalet
and Bengio, 2005) that controls the kurtosis of the distribution. In this
work, we consider using an annealed $\tau$, similar to Equation 2 where
$\tau\to 2.5$ and starts with $\tau=0.01$. This is done to allow the model to
avoid large gradient variance in updates early on. In the context of using the
Gumbel-Softmax in LVMS, this allows the model to become robust to non-greedy
actions gradually throughout training, we would expect such exploration to
improve generalization proportional to the vocabulary size.

##  4 Experimental Setup

We carry out experiments for a 2-hidden layer Long-Short Term Memory (LSTM)
model with embedding size $|e|=400$, Backpropogation Through Time (BPTT)
length $35$ and variational dropout Gal and Ghahramani (2016) with rate
$p_{d}$=0.2 for input, hidden and output layers. The ECOC-NSP model is trained
using Binary Cross Entropy loss as shown in Equation 6, where $k$ is a group
error-checking codewords corresponding to a codeword $C$. The gradients can
then expressed as
$\frac{\delta\mathcal{L}}{\delta\theta}=(y−\sigma(h^{T}\theta))\cdot h^{T}$.

$\begin{gathered}\mathcal{L}_{\theta}=\max_{k}\prod_{c}^{C}\big{[}y_{c}\log\big{(}\sigma_{c}(h^{T}\theta\big{)}+\\\
(1-y_{c})\log\big{(}1-\sigma_{c}(h^{T}\theta)\big{)}\big{]}\end{gathered}$ (6)

Baselines for ECOC-NSP The first set of experiments include comparisons
against the most related baselines, which include Sample-Softmax (Bengio et
al., 2003; Bengio and Senécal, 2008), Hierarchical Softmax (HS), AS (Grave et
al., 2016), and NCE (Mnih and Teh, 2012). For HS, we use a 2-hidden layer tree
with a branching factor (number of classes) of $\sqrt{|\mathcal{V}|}$ by
default. For AS, we split the output into 4 groups via the unigram
distribution (in percentages of total words 5%-15%-30%-100%). For NCE, we set
the noise ratio to be 0.1 for PTB and 0.2 for WikiText-2 and WikiText-103.
Training is carried out until near convergence ($\epsilon\approx 40$), the
randomly initialized HS and Sampled Softmax of which take longer
($\epsilon\in$ [55-80]). Table 1 reports the results for
$\log_{2}|\mathcal{V}|^{2}$ number of samples in the case of Rand/Uni-Sample-
SM. For Rand/Unigram Hierarchical SM, we use a 2-hidden layer tree with 10
classes per child node.

Baselines for ECOC-NSP Mixture Sampling To test Latent Variable Mixture
Sampling (LVMS), we directly compare its application in HS and ECOC, two
closely related latent NSP methods. Additionally, we compare the performance
of LVMS against the most related sampling-based supervised learning technique
called scheduled sampling (SS) (Bengio et al., 2015). For SS with cross-
entropy based training (SS-CE), we also consider using a baseline of the soft-
argmax (Soft-SS-CE) where a weighted average embedding is generated
proportional to the predicted probability distribution.

Evaluation Details To compute perplexities for ECOC-NSP, the codewords are
viewed in terms of a product of marginal probabilities. Hence, when
approximating the posterior, we treat latent codes as a factorial distribution
Marlow et al. (1965) where each bit is independent of one another ($c_{i}\perp
c_{j}\;\forall i,j$) . Among the $k_{i}$ error checks corresponding to a
particular tokens codeword $C_{k_{i}}$, we choose the most probable of these
checks, as shown in Equation 7, when computing the binary cross-entropy loss
at training time.

$p_{\hat{C}}=\max_{k_{i}}\prod_{c_{k_{i}}}^{C_{k_{i}}}\Big{(}\phi(x_{t},h_{t})\Big{)},\;i=[0,1..,k]$
(7)

At test time, if the predicted codeword $\hat{C}$ falls within the $k$ error-
checking bits of codeword $C$, it is deemed a correct prediction and assigned
the highest probability of all $k$ predictions. Note that we only convert the
ECOC predictions to perplexities to be comparable against baselines (we could
use Hamming Distance or Mean Reciprocal Rank when the codes are easily decoded
by ordered semantically or by Hamming distance).

Figure 2: ECOC-NSP Peplexity vs. Decoder Parameters (corresponding to 14/20/40
codeword bits for Penn-TreeBank and 17/40/100 codeword bits for
WikiText-2/103)

##  5 Results

##### Error-Correcting Output Coded NSP

We first compare our proposed ECOC-NSP to aforementioned methods that
approximate softmax normalization, using binary trees and latent codes that
are ordered according to unigram frequency (Uni-Hierarchical-SM and Uni-ECOC).
This is also the same ordering we use to compare our proposed CMS-ECOC
sampling method to scheduled sampling (Bengio et al., 2015) in standard cross-
entropy training with softmax normalization. Although, these are not directly
comaprable, since ECOC-NSP introduces a whole new paradigm, we use the common
evaluation measures of Hamming distance and accuracy to have some kind of
baseline with which we can compare our proposed method to. Equation 2 shows
how the reduction in perplexity as the number of ECOC-LSTM decoder parameters
increase as more bits are added to the codeword. For PTB, large perplexity
reductions are made between 14-100 codebits, while between 100-1000 codebits
there is a gradual decrease. In contrast, we see that there is more gained
from increasing codeword size for WikiText-2 and WikiText-103 (which preserve
the words that fall within the long-tail of the unigram distribution). We find
the discrepancy in performance between randomly assigned codebooks and ordered
codebooks is more apparent for large compression ($|C|<|\mathcal{V}|/10$).
Intuitively, the general problem of well-separated codes is alleviated as more
bits are added.

Equation 1 shows that overall ECOC with a rank ordered embedding similarity
$\mathcal{C}$ (Embedding-ECOC) almost performs as well as the full-softmax
(8.02M parameters) while only using 1000 bits for PTB ($|\mathcal{V}|/20$ and
) and 5K bits for WikiText-2 ($|\mathcal{V}|/25$) and WikiText-103
($|\mathcal{V}|/30$). The HS-based models use a 2-hidden layer tree with 10
tokens per class, resulting in 4.4M parameters for PTB, 22.05M parameters for
WikiText-2 (full softmax - 40.1M) and WikiText-103. Moreover, we find there is
a consistent improvement in using Embedding-ECOC over using a random codebook
(Random-ECOC) and a slight improvement over using a unigram ordered codebook
(Uni-ECOC). Note that in both Embedding-ECOC and Uni-ECOC, the number of
error-checking bits are assigned inversely proportional to the rank position
when ordering embedding similarities (as discussed in Equation 3.1.1) and
unigram frequency respectively. We also found that too many bits e.g
$|C|=|\mathcal{V}|$ takes much longer ($\epsilon\in$ [20-30] more for PTB) to
converge with negligible perplexity reductions. Hence, the main advantage of
ECOC-NLVMS is the large compression rate while maintaining performance (e.g
PTB with $|C|=40$, there is less than 2 perplexity points compared to the full
softmax).

Model | PTB | WikiText-2 | WikiText-103
---|---|---|---
| Val. | Test | Val. | Test | Val. | Test
Full SM | _86.19_ | _79.24_ | _124.01_ | _119.30_ | _56.72_ | _49.35_
Rand-Sample-SM | 92.14 | 81.82 | 136.47 | 129.29 | 68.95 | 59.34
Uni-Sample-SM | 90.37 | 81.36 | 133.08 | 127.19 | 66.23 | 57.09
Rand-Hierarchical-SM | 94.31 | 88.50 | 133.69 | 127.12 | 62.29 | 54.28
Uni-Hierarchical-SM | 92.38 | 86.70 | 130.26 | 124.83 | 62.02 | 54.11
Adaptive-SM | 91.38 | 85.29 | 118.89 | 120.92 | 60.27 | 52.63
NCE | 96.79 | 89.30 | 131.20 | 126.82 | 61.11 | 54.52
Random-ECOC | 91.00 | 87.19 | 131.01 | 123.29 | 56.12 | 52.43
Uni-ECOC | 86.44 | 82.29 | 129.76 | 120.51 | _52.71_ | _48.37_
Embedding-ECOC | _84.40_ | _77.53_ | _125.06_ | _120.34_ | 57.37 | 49.09

Table 1: LSTM Perplexities for Full Softmax (SM), Sample-Based SM (Sample-SM),
Hierarchical-SM (HSM), Adaptive-SM, NCE and ECOC-NSP.

##### Latent Variable Mixture Sampling Text Generation Results

| B1 | B2 | B3 | B4 | R-L | MET
---|---|---|---|---|---|---
Full-SM | 71.09 | 51.33 | 32.85 | 24.67 | 50.28 | 52.70
SS-SM | 73.23 | 52.81 | 33.37 | 26.11 | 52.60 | 54.51
Soft-SS-SM | _73.54_ | _53.01_ | _33.26_ | _27.13_ | _54.49_ | _54.83_
SS-Adaptive-SM | 70.45 | 50.22 | 31.38 | 23.59 | 51.88 | 51.83
SS-Hierarchical-SM | 67.89 | 48.42 | 30.37 | 22.91 | 49.39 | 50.48
CLVMS-Hierarchical-SM | 69.70 | 49.52 | 31.91 | 24.19 | 51.35 | 51.20
DLVMS-Hierarchical-SM | 71.04 | 50.61 | 32.26 | 24.72 | 52.83 | 52.36
SS-ECOC | 72.02 | 52.03 | 32.57 | 25.42 | 51.39 | 53.51
Soft-SS-ECOC | 72.78 | 53.29 | 33.15 | 25.93 | 52.07 | 54.22
CLVMS-ECOC | _74.70_ | _53.09_ | _34.28_ | _27.05_ | _53.67_ | _55.62_
DLVMS-ECOC | _74.92_ | _53.56_ | _34.70_ | _27.81_ | _54.02_ | _55.85_

Table 2: MSCOCO Test Results on BLEU (B), ROUGE-L (R-L) & METEOR (MET) Eval.
Metrics

Table 2 shows all results of LVMS when used in HS and ECOC-based NSP models
for the MSCOCO image captioning dataset Lin et al. (2014) using the Karpathy
validation and test splits using a beam search with a width of 5 at test time.
We use $|c|=200$ to account for vocabulary size $|V|=10^{3}$, leaving
$|c|-\log_{2}(|V|)=186$ error-check bits leftover $\forall C\in\mathcal{C}$.
The HS uses the Categorical Concrete distribution for DLVMS-HS and Binary
Concrete Distribution for DCMS-ECOC. In our experiments we found $\tau=2.0$ to
be the upper threshold from an initial grid search of
$\tau\in\\{0.2,0.5,1,2,5,10\\}$ for both DLVMS-HS and DCMS-ECOC, where
$\tau<2$ corresponds to little exploration and $\tau>2$ results in too much
exploration, particularly early on when the model is performing larger
gradient updates. CLVMS-Hierarchical-SM and CLVMS-ECOC both monotonically
increases $\tau$ according to Equation 2.

Both HS and ECOC use Embedding ordered decoder matrix (we omit the -Embedding
extension). We baseline this against both SS and the soft-argmax version of
SS, the most related sample-based supervised learning approach to LVMS.
Furthermore, we report results on CLVMS-ECOC (Curriculum-LVMS ECOC) which
mixes between true targets and codewords predictions according to the schedule
in Equation 2 and a differentiable extension of LVMS via samples from the
Gumbel-Softmax (DCMS-ECOC). For both DCMS-ECOC and DLVMS-Hierarchical-SM we
sample from each softmax defined along the path to the target code at training
time. We find that using a curriculum in CLVMS-ECOC with a semantically
ordered codebook outperforms the full softmax that uses scheduled sampling
(SS-SM) and its weighted-variant (Soft-SS-SM). Moreover, DLVMS-ECOC further
improves over CLVMS-ECOC on MSCOCO. LVMS makes a consistent improvement over
SS. This suggests that LVMS is an effective alternative for latent variable
based NSPs in particular (mixture sampling is ill-suited to one-hot targets as
they are extremely sparse).

##  6 Conclusion

This work proposed an error-correcting neural language model and a novel
Latent Variable Mixture Sampling method for latent variable models. We find
that performance is maintained compared to using the full conditional and
related approximate methods, given a sufficient codeword size to account for
correlations among classes. This corresponds to 40 bits for PTB and 100 bits
for WikiText-2 and WikiText-103. Furthermore, we find that performance is
improved when rank ordering the codebook via embedding similarity where the
query is the embedding of the most frequent word.
Lastly, we introduced Latent Variable Mixture Sampling to mitigate exposure
bias. This can be easily integrated into training latent variable-based
language models, such as the ECOC-based language model. We show that this
method outperforms the well-known scheduled sampling method with a full
softmax, hierarchical softmax and an adaptive softmax on an image captioning
task, with less decoder parameters than the full softmax with only 200 bits,
2% of the original number of output dimensions.

## References

  * Bengio et al. (2015) Samy Bengio, Oriol Vinyals, Navdeep Jaitly, and Noam Shazeer. 2015.  Scheduled sampling for sequence prediction with recurrent neural networks.  In _Advances in Neural Information Processing Systems_ , pages 1171–1179.
  * Bengio and Senécal (2008) Yoshua Bengio and Jean-Sébastien Senécal. 2008.  Adaptive importance sampling to accelerate training of a neural probabilistic language model.  _IEEE Transactions on Neural Networks_ , 19(4):713–722.
  * Bengio et al. (2003) Yoshua Bengio, Jean-Sébastien Senécal, et al. 2003.  Quick training of probabilistic neural nets by importance sampling.  In _AISTATS_ , pages 1–9.
  * Berger (1999) Adam Berger. 1999.  Error-correcting output coding for text classification.  In _IJCAI-99: Workshop on machine learning for information filtering_.
  * Chen et al. (2015) Welin Chen, David Grangier, and Michael Auli. 2015.  Strategies for training large vocabulary neural language models.  _arXiv preprint arXiv:1512.04906_.
  * Gal and Ghahramani (2016) Yarin Gal and Zoubin Ghahramani. 2016.  A theoretically grounded application of dropout in recurrent neural networks.  In _Advances in neural information processing systems_ , pages 1019–1027.
  * Goodman (2001) Joshua T Goodman. 2001.  A bit of progress in language modeling.  _Computer Speech & Language_, 15(4):403–434.
  * Goyal et al. (2017) Kartik Goyal, Chris Dyer, and Taylor Berg-Kirkpatrick. 2017.  Differentiable scheduled sampling for credit assignment.  _arXiv preprint arXiv:1704.06970_.
  * Grandvalet and Bengio (2005) Yves Grandvalet and Yoshua Bengio. 2005.  Semi-supervised learning by entropy minimization.  In _Advances in neural information processing systems_ , pages 529–536.
  * Grave et al. (2016) Edouard Grave, Armand Joulin, Moustapha Cissé, David Grangier, and Hervé Jégou. 2016.  Efficient softmax approximation for gpus.  _arXiv preprint arXiv:1609.04309_.
  * Hamming (1950) Richard W Hamming. 1950.  Error detecting and error correcting codes.  _Bell System technical journal_ , 29(2):147–160.
  * Jang et al. (2016) Eric Jang, Shixiang Gu, and Ben Poole. 2016.  Categorical reparameterization with gumbel-softmax.  _arXiv preprint arXiv:1611.01144_.
  * Lin et al. (2014) Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick. 2014.  Microsoft coco: Common objects in context.  In _European conference on computer vision_ , pages 740–755. Springer.
  * Maddison et al. (2016) Chris J Maddison, Andriy Mnih, and Yee Whye Teh. 2016.  The concrete distribution: A continuous relaxation of discrete random variables.  _arXiv preprint arXiv:1611.00712_.
  * Marlow et al. (1965) William H Marlow et al. 1965.  Factorial distributions.  _The Annals of Mathematical Statistics_ , 36(3):1066–1068.
  * Mnih and Teh (2012) Andriy Mnih and Yee Whye Teh. 2012.  A fast and simple algorithm for training neural probabilistic language models.  _arXiv preprint arXiv:1206.6426_.
  * Morin and Bengio (2005) Frederic Morin and Yoshua Bengio. 2005.  Hierarchical probabilistic neural network language model.  In _Aistats_ , volume 5, pages 246–252.
  * Oda et al. (2017) Yusuke Oda, Philip Arthur, Graham Neubig, Koichiro Yoshino, and Satoshi Nakamura. 2017.  Neural machine translation via binary code prediction.  _arXiv preprint arXiv:1704.06918_.
  * Sejnowski and Rosenberg (1987) Terrence J Sejnowski and Charles R Rosenberg. 1987.  Parallel networks that learn to pronounce english text.  _Complex systems_ , 1(1):145–168.
  * Shi and Yu (2018) Kaiyu Shi and Kai Yu. 2018.  Structured word embedding for low memory neural network language model.  _Proc. Interspeech 2018_ , pages 1254–1258.
  * Shu and Nakayama (2017) Raphael Shu and Hideki Nakayama. 2017.  Compressing word embeddings via deep compositional code learning.  _arXiv preprint arXiv:1711.01068_.
  * Sutskever et al. (2014) Ilya Sutskever, Oriol Vinyals, and Quoc V Le. 2014.  Sequence to sequence learning with neural networks.  In _Advances in neural information processing systems_ , pages 3104–3112.

——


Corresponding author<EMAIL_ADDRESS>
#  Universal finite-time scaling in the transcritical, saddle-node, and
pitchfork discrete and continuous bifurcations

Álvaro Corral  Centre de Recerca Matemàtica, Edifici C, Campus Bellaterra,
E-08193 Barcelona, Spain  Departament de Matemàtiques, Facultat de Ciències,
Universitat Autònoma de Barcelona, E-08193 Barcelona, Spain  Complexity
Science Hub Vienna, Josefstädter Stra$\beta$e 39, 1080 Vienna, Austria

###### Abstract

Bifurcations are one of the most remarkable features of dynamical systems.
Corral et al. [Sci. Rep. 8(11783), 2018] showed the existence of scaling laws
describing the transient (finite-time) dynamics in discrete dynamical systems
close to a bifurcation point, following an approach that was valid for the
transcritical as well as for the saddle-node bifurcations. We reformulate
those previous results and extend them to other discrete and continuous
bifurcations, remarkably the pitchfork bifurcation. In contrast to the
previous work, we obtain a finite-time bifurcation diagram directly from the
scaling law, without a necessary knowledge of the stable fixed point. The
derived scaling laws provide a very good and universal description of the
transient behavior of the systems for long times and close to the bifurcation
points.

##  I Introduction

Bifurcations involve a sudden and qualitative change in the asymptotic
solutions of a dynamical system when one parameter (or more) is gradually
modified Strogatz_book . This implies that bifurcations are able to explain
profound variations in the stationary state of many real systems when some
conditions (given by the bifurcation parameter) are slightly different. These
real systems include physical, chemical, biological, environmental, and any
sort of system whose evolution can be formulated in terms of ordinary
differential equations (or their discrete counterpart). Remarkable examples
are lasers, oscillating chemical reactions, activity of genes, animal or plant
populations, and epidemic spreading. Strogatz_book ; Sole_phase_transitions .

Trying to strength the connection between branching processes Harris_original
; branching_biology ; Corral_FontClos  and thermodynamic phase transitions
Stanley ; Yeomans1992 ; Stanley_rmp , as well as to better understand the
emergence of critical behavior in the infinite system-size limit, Garcia-
Millan et al. GarciaMillan ; Corral_garciamillan  studied the renowned Galton-
Watson process. They analytically established the existence of finite-size
scaling Privman  and universality Stanley_rmp  in the relation between the
probability of survival at different finite times and the mean offspring
number per individual when the offspring distribution is arbitrary (with the
only requirement of a finite variance; the opposite case has been studied more
recently and leads to a continuum of universality classes Falco_Corral ).
Thus, the Galton-Watson process displays a sudden change analogous to a phase
transition. Notably, the aforementioned authors GarciaMillan ;
Corral_garciamillan  found an exact analytical expression for the scaling
function appearing in the scaling law describing the probability of survival.
Other aspects of the finite-size properties of the Galton-Watson process were
analyzed in Ref. Corral_garcia_moloney_font .

The same authors GarciaMillan ; Corral_garciamillan , however, did not pay
attention to the fact that there is a transcritical bifurcation underlying the
phase transition in the Galton-Watson process, as the time evolution of the
probability of survival is governed by a discrete dynamical system. Reference
Corral_Alseda_Sardanyes  put the results of Garcia-Millan et al. GarciaMillan
; Corral_garciamillan  in the context of bifurcations in dynamical systems,
showing that the results were valid beyond the Galton-Watson process, as they
arose from a universal, previously unnoticed scaling characteristic of the
transcritical bifurcation. It was also shown Corral_Alseda_Sardanyes  that
this property was present in the saddle-node bifurcation as well. It is
noteworthy how, in the approach of Ref. Corral_Alseda_Sardanyes , finite-size
scaling in the Galton-Watson process translates into finite-time scaling in
the framework of dynamical systems (understanding the independent variable of
the system as time, yielding a description of the transient dynamics), due to
the fact that in the former process, system size (given by the depth of the
underlying branching tree) and time (given by the number of generations) are
so intertwined that they are in fact the same thing (although the evolution of
dynamical systems is more naturally described in terms of time).

In the present article we are guided by this parallelism between phase
transitions and bifurcations in dynamical systems, summarized in Table 1. We
provide a clear case in which statistical-physics concepts can be used to
better understand dynamical systems, concretely, bifurcations. Despite of the
fact that our approach looks similar to that in Ref. Corral_Alseda_Sardanyes ,
we expand the results of that work in several directions. In Sec. 2 we briefly
summarize the results of Ref. Corral_Alseda_Sardanyes . Section 3 deals with
discrete dynamical systems, for which we derive and test a scaling law for the
distance of the solution to a fixed point that can be stable or unstable, and
when this is constant (independent on the bifurcation parameter), we directly
obtain a “finite-time bifurcation diagram” that is valid for a variety of
bifurcations, including the pitchfork one (in addition to the transcritical
and saddle-node cases studied in Ref. Corral_Alseda_Sardanyes ). In Sec. 4 we
extend the previous results Corral_Alseda_Sardanyes  as well as our new
results to continuous dynamical systems, whereas in the subsequent section we
study an alternative scaling law for the distance of the solution to a stable
fixed point, exploring also the limits of validity of these scaling laws.
Finally, after discussing the results, an Appendix shows how the well-known
phenomenon of critical slowing down is contained in our scaling laws, allowing
the derivation of the exponents of the temporal power-law decay to the
attractor at the bifurcation point as well as the values of the characteristic
decay times outside this point.

Table 1:  Parallelisms between thermodynamic phase transitions, branching processes, and bifurcations in dynamical systems (thermodynamic phase transitions are illustrated with the example of a magnetic material at zero external field). The notation we use in the main text is introduced in the third column. In dynamical systems, fixed points are sometimes referred to as critical points, but we do not use such terminology as there is no correspondence with the usage of the same term in phase transitions.  Phase transitions | Branching processes | Dynamical systems
---|---|---
Control parameter | Mean offspring number | Bifurcation parameter
(as temperature $T$) | $\langle N\rangle$ | $\mu$
Order parameter | Probability of survival (at infinity) | Fixed-point attractor
in thermodynamic limit | $\mathcal{P}$ | (or stable fixed point)
(as magnetization $\mathbb{M}$) |  | $q$
Critical point | Critical point | Bifurcation point
(as $(T^{*},\mathbb{M}^{*})$ with $\mathbb{M}^{*}=0$) |  $(\langle N\rangle^{*},\mathcal{P}^{*})$ with $\langle N\rangle^{*}=1$, $\mathcal{P}^{*}=0$ |  $(\mu^{*}$, $q^{*})$
System size | Number of generations | Independent variable
|  | (time $n$, $\ell$, $t$, $\tau$)
Finite-size scaling | Finite-size scaling $=$ | Finite-time scaling
| finite-time scaling |

##  II Previous research: The discrete transcritical and saddle-node
bifurcations

Let us consider the following one-dimensional discrete dynamical system or map
Strogatz_book

$x_{n+1}=f(x_{n})$ (1)

where $n=0,1\dots$ denotes a discrete time (counting time generations) and
$x_{n}$ are real numbers; thus, $f(x)$ (also called a map) is a real function.
It is assumed that the system has a linearly stable fixed point (also referred
to here as an attractive fixed point) at $x=q$, with $q=f(q)$ and
$|f^{\prime}(q)|<1$ (the prime denotes differentiation and the bars absolute
value). It is implicit that the function $f(x)$ depends on a bifurcation
parameter $\mu$ and so, $q$ may depend on $\mu$ as well.

Reference Corral_Alseda_Sardanyes  (based on results from Ref. GarciaMillan )
showed that in such a discrete dynamical system, for large values $\ell$ of
the discrete time, close to a bifurcation point, and when the initial
condition $x_{0}$ belongs to the basin of attraction of the fixed point, the
difference between the stable fixed point $q$ and the solution at time $\ell$
fulfills a finite-time scaling law,

$q-x_{\ell}=\frac{1}{D^{*}\ell}\,G(z),$ (2)

where the rescaled variable $z$ is defined as

$z=\ell(\tilde{m}-1),$ (3)

with $\tilde{m}=f^{\prime}(q)$, and $D^{*}=f^{\prime\prime}(q^{*})/2\neq 0$,
with $q^{*}$ denoting the value that the fixed point takes at the bifurcation
point. Remarkably, the expression does not depend on the initial condition.
Even more remarkably, the scaling function $G(z)$ turned out to be independent
on $f(x)$ and therefore “universal,” being

$G(z)=\frac{ze^{z}}{e^{z}-1}.$ (4)

It is worth mentioning that Eq. (2), together with Eq. (3) constitutes a clear
example of what is called a scaling law in statistical physics
Christensen_Moloney ; Corral_csf . Other authors use the term scaling law to
refer to (univariate) power-law functions; certainly, power laws are a
particular case of scaling laws (when the scaling function is constant);
nevertheless, in this article a scaling law refers to something akin to Eqs.
(2) plus (3) and not to a simple power law.

Note that, as the fixed point $q$ is linearly stable, then $\tilde{m}<1$, and
$z<0$. Restricted to $z\leq 0$ (extending $z$ to the bifurcation point $z=0$,
where the attractor loses its linear stability) the function $G(z)$ is
monotonically increasing with $z$, with a maximum at $z=0$. This indicates
that the closer to the bifurcation point, the longer it takes to reach a small
neighborhood of the attractor; this is quantified by the scaling law, and the
rescaled variable $z$ provides the appropriate way to compute the “distance”
to the bifurcation point (depending on the value of $\ell$).

In consequence, the attractors of certain one-dimensional discrete dynamical
systems given by Eq. (1) are approached in the same universal way, with the
only condition that $f^{\prime\prime}(q^{*})\neq 0$. This was found to be
valid for the transcritical and the saddle-node bifurcations
Corral_Alseda_Sardanyes . Knowing the dependence of $q$ on the bifurcation
parameter, it was also possible to find scaling laws for the relation between
$x_{\ell}$ and $\tilde{m}$ (eliminating $q$ from Eq. (2)) in the transcritical
and saddle-node bifurcations (for large $\ell$ and close to the bifurcation
point). Nevertheless, this second scaling laws were not fully universal, in
the sense that they were different for each bifurcation (although the scaling
functions for each one were mathematically related through simple
transformations).

##  III  Finite-time bifurcation diagram in discrete bifurcations

###  III.1 Theoretical derivation

We start with the one-dimensional discrete dynamical system considered in Ref.
Corral_Alseda_Sardanyes  and given by Eq. (1). An important difference between
our approach and that previous work is that, for us, the reference fixed point
can be either linearly stable or linearly unstable (or neutral), and then,
instead of by $q$, it will be denoted by $p$, verifying $|f^{\prime}(p)|<1$ or
$|f^{\prime}(p)|>1$ (or $|f^{\prime}(p)|=1$), respectively. In other words,
the stability of the fixed point does not matter for us. Also, the previous
work Corral_Alseda_Sardanyes  assumed, as mentioned,
$D^{*}=f^{\prime\prime}(q^{*})/2\neq 0$, which holds for the transcritical and
saddle-node bifurcations, but not in other cases.

We assume here that $f^{\prime\prime}(p^{*})$ and $f^{\prime\prime}(p)$ may be
zero and define $P=-f^{(k)}(p)/k!$, with $k\geq 2$ the lowest order of the
derivative for which $f^{(k)}(p)\neq 0$, and with $k$ in parenthesis denoting
the $k-$th order derivative (with $k$ integer; the case studied in Ref.
Corral_Alseda_Sardanyes  corresponds to $k=2$). Given the initial condition
$x_{0}$, the $n-$th iteration of $f$ is $f(\dots
f(f(x_{0}))\dots)=f^{n}(x_{0})$. For sufficiently large $n$, we perform a
Taylor expansion of $f(f^{n}(x_{0}))$ around the (stable or unstable) fixed
point $p$, so

$f^{n+1}(x_{0})=p+m(f^{n}(x_{0})-p)-P(f^{n}(x_{0})-p)^{k}+\mathcal{O}(f^{n}(x_{0})-p)^{k+1},$
(5)

with $m=f^{\prime}(p)$. If $p$ is attractive, as $n$ is large, it is obvious
that $f^{n}(x_{0})$ will be close to $p$; however, when $p$ is unstable (that
is, repulsive), we assume that there will be an attractive fixed point $q$ in
between $x_{0}$ and $p$, and $f^{n}(x_{0})$ will be close to $q$; moreover,
close to the bifurcation point, $q$ will be close to $p$ and therefore
$f^{n}(x_{0})$ will be close to $p$. It is implicit then that the bifurcation
is continuous. It is assumed also that $x_{0}$ belongs to the basins of
attraction of the attractive fixed point.

Introducing the “distance” (which can be negative) between the solution at
time $n$ and the fixed point $p$ as $d_{n}=f^{n}(x_{0})-p$ we can write Eq.
(5) as $d_{n+1}=md_{n}-Pd_{n}^{k}$, up to order $k$. It is convenient to
define $c_{n}=1/d_{n}^{\alpha}$, with $\alpha$ arbitrary footnote_aina , and
from here, as $d_{n}$ will be small,

$c_{n+1}=\frac{c_{n}}{m^{\alpha}}+\frac{\alpha
Pd_{n}^{k-1-\alpha}}{m^{\alpha+1}}$

plus higher-order terms Aina_report . The expression is particularly simple
taking $\alpha=k-1$, then, iterating $\ell$ times we arrive at

$c_{n+\ell}\simeq\frac{c_{n}}{m^{(k-1)\ell}}+\frac{(k-1)P}{m^{k}}\left(1+\frac{1}{m^{k-1}}+\dots+\frac{1}{m^{(k-1)(\ell-1)}}\right)$
$=\frac{c_{n}}{m^{(k-1)\ell}}+\frac{(k-1)P}{m^{k}}\frac{m^{(k-1)\ell}-1}{m^{(k-1)(\ell-1)}(m^{k-1}-1)}.$

We deal with the existence of a bifurcation at $f^{\prime}(p)=1$ (ignoring
bifurcations given by $f^{\prime}(p)=-1$); then, similarly as in Eq. (3),
Corral_Alseda_Sardanyes , we define,

$y=\ell(m-1),$ (6)

with $m=f^{\prime}(p)$ and therefore $m^{\ell}\rightarrow e^{y}$ in the limit
of large $\ell$. From here,

$c_{n+\ell}\simeq\frac{c_{n}}{e^{(k-1)y}}+(k-1)P\,\frac{e^{(k-1)y}-1}{e^{(k-1)y}(k-1)y/\ell}.$

We may define $\ell^{\prime}=n+\ell$ but for $\ell\gg n$ we have
$\ell\simeq\ell^{\prime}$. In that case, the second term dominates (it is
proportional to $\ell$), and

$c_{\ell}\simeq\frac{(k-1)P^{*}\ell}{G((k-1)y)},$

with $G(y)$ given by Eq. (4) and approximating (to the lowest order) $P\simeq
P^{*}$ (the value at the bifurcation point, as for large $\ell$ and finite
$y$, the parameter $m$ has to be close to one and the system has to be close
to the bifurcation).

Coming back to the distance between the solution and the fixed point, we have
that $d_{\ell}^{k-1}=1/c_{\ell}$ and then

$\left(f^{\ell}(x_{0})-p\right)^{k-1}={\frac{1}{(k-1)P^{*}\ell}\,G((k-1)y)},$
(7)

which constitutes, together with Eq. (6), a finite-time scaling law. From now
on we write it as an equality, although it is only exact in the limits
$\ell\rightarrow\infty$ and $m\rightarrow 1$. This scaling law contains the
one in Ref. Corral_Alseda_Sardanyes  as a special case, taking $k=2$ and
considering that the fixed point is stable (changing $p$ for $q$). It is
remarkable that we obtain the same scaling function $G(y)$ as in Refs.
GarciaMillan ; Corral_Alseda_Sardanyes , despite our problem is more general.
Note that we use $y$ instead of $z$, as $z$ in Ref. Corral_Alseda_Sardanyes
was restricted to $z\leq 0$, but as $p$ can now be unstable we do not have
such a restriction. Moreover, in our case, $m$ can be considered to play the
role of an alternative, “natural” bifurcation parameter (and therefore has to
be related to $\mu$) as we do not have the restriction $\tilde{m}<1$. This
means that $m$ quantifies the distance to the bifurcation point (given by
$m^{*}=1$). We also realize that, when $k$ is odd, $P^{*}$ has to be larger
than zero (for the scaling law to hold), as the scaling function $G$ is
positive for any finite $y$ but $x_{\ell}-p$ can take both positive and
negative values. In contrast, when $k$ is even there is no restriction on the
sign of $P^{*}$, and $x_{\ell}-p$ will have the same sign as $P^{*}$;
therefore, in this case the scaling law is only valid at one side in the
bifurcation diagram. It is straightforward but very useful to realize that the
scaling law (7) can also be written as

$f^{\ell}(x_{0})-p=(\pm
1)^{k}\left({\frac{1}{(k-1)P^{*}\ell}\,G((k-1)y)}\right)^{1/(k-1)},$ (8)

where the term $(\pm 1)^{k}$ arises from the real values of $\sqrt[k-1]{1}$.

###  III.2 Bifurcation arising from the smoothness of the scaling function

Using the obtained scaling law, we can see how the sharp change represented by
the bifurcation emerges in the limit $\ell\rightarrow\infty$ from the
asymptotic properties of a continuous and smooth scaling function $G(y)$.
Indeed, from Eq. (4),

$G(y)\rightarrow\left\\{\begin{array}[]{ll}0&\mbox{ for
}y\rightarrow-\infty,\\\ 1&\mbox{ for }y\rightarrow 0,\\\ y&\mbox{ for
}y\rightarrow\infty,\\\ \end{array}\right.$

and substituting in the scaling law, Eq. (8),

$f^{\ell}(x_{0})-p\rightarrow\left\\{\begin{array}[]{lll}0&\mbox{ for
}m<1&\mbox{ and }\ell\rightarrow\infty,\\\ (\pm
1)^{k}\sqrt[k-1]{\frac{1}{(k-1)P^{*}\ell}}&\mbox{ for }m=1,&\\\ (\pm
1)^{k}\sqrt[k-1]{\frac{m-1}{P^{*}}}&\mbox{ for }m>1&\mbox{ and
}\ell\rightarrow\infty,\\\ \end{array}\right.$ (9)

with the obvious meaning that when $m<1$, the fixed point $p$ is (linearly)
stable and the solution tends to it for $\ell\rightarrow\infty$, whereas when
$m>1$, the fixed point $p$ is unstable and the solution tends to another fixed
point $q$.

In the latter case, the value of the stable fixed point $q$ (minus $p$) is
obtained from the Taylor expansion around the unstable fixed point $p$,
without any knowledge of $q$. A case of great interest is when $p$ is zero, or
constant (for all the considered values of $\mu$), for which the scaling law
yields the equivalent of the bifurcation diagram but for finite times, without
the need of calculating the stable fixed point $q$; i.e., only the knowledge
of the constant (stable or unstable) fixed point $p$ is required.

In the previous subsection we have seen that, when $k$ is odd, then $P^{*}>0$
and $x_{\ell}-p$ can take any sign. Now we can see, from Eq. (9), which yields
$q-p$ (with $q=p$ for $m<1$), that $x_{\ell}>q$ when $q>p$ and $x_{\ell}<q$
when $q<p$. The reason is that $G(y)>y$ when $y>0$ and $G(y)>0$ when $y<0$ (in
fact, always). Thus, given an attractive fixed point $q$, the scaling law is
only valid for initial conditions at one of its sides. The same is true when
$k$ is even (but this is already clear from the previous subsection).

###  III.3 Numeric verification

For a given dynamical system, characterized by a value of $k$, the fulfillment
of the scaling law (8) for any value of the discrete time $\ell$ (large) and
$m$ (close to one, related to the bifurcation parameter) must become apparent
displaying $\sqrt[k-1]{(k-1)P^{*}\ell\,}\,(x_{\ell}-p)$ versus $y=\ell(m-1)$,
yielding the shape of the function $(\pm 1)^{k}[G((k-1)y)]^{1/(k-1)}$, which
does not depend neither on $\ell$ nor on $m$ individually, but on
$\ell(m-1)=y$ (note the metric factor $k-1$ multiplying $y$ in $G()$ Privman
). We obtain in this way a law of corresponding states Stanley  that signals
the existence of universality, as only $k$ and the first non-zero derivatives
$m$ and $P^{*}$ of the map $f(x)$ matter. However, there is an additional
level of universality, as the same scaling law is valid for different
bifurcations, taking into account that the value of $k$ may change from one
bifurcation to another. In that case, for different $k$, Eq. (7) indicates
that one should display the more general rescaling
${(k-1)P^{*}\ell\,(x_{\ell}-p)^{k-1}}$ versus $(k-1)\ell(m-1)$, leading to a
collapse of all the data into the unique function $G$.

We can test the validity of the scaling law with specific examples. Figures 1
and 2 analyze, close to the bifurcation point, the behavior of several
discrete dynamical systems given by

$f(x)=(1+\mu)x-x^{k},$ (10)

where the bifurcation parameter is $\mu$ and $k=2,3...$ corresponds to
different bifurcations. In this way, $k=2$ yields the normal form of the
transcritical bifurcation and $k=3$ gives the normal form of the (so-called
supercritical) pitchfork bifurcation Strogatz_book . Independently of the
value of $k$, all the dynamical systems represented by this $f(x)$ have $p=0$
as a fixed point (it does not matter if stable or unstable) for all $\mu$.
Then, $m=f^{\prime}(p)=1+\mu$ (thus, all the bifurcations for different values
of $k$ take place at $\mu=0$), $y=\ell(m-1)=\mu\ell$ and $P=-f^{(k)}(p)/k!=1$.
Figures 1(a) and 2(a) display, for $k=3$ and $4$ respectively, the solution
$f^{\ell}(x_{0})$ as a function of the bifurcation parameter $\mu$ for
different values of $\ell$ (and changing also $x_{0}$). What we obtain are the
usual profiles that approach the underlying bifurcation diagram in the limit
$\ell\rightarrow\infty$.

Figures 1(b) and 2(b) display essentially the same data but applying the
rescaling given by the scaling law, under the form directly suggested by Eq.
(8). The observed nearly-perfect data collapses clearly support the validity
of the scaling law. Figure 2(c) gives further support, including, in addition,
results from $k=2$ to $k=5$, as well as from the logistic map at its
transcritical bifurcation, using the more general rescaling arising directly
from Eq. (7). This rescaling, with the “generalized distance”
$(x_{\ell}-p)^{k-1}$ to the fixed point, leads to a very good data collapse;
this indicates a strong universal behavior for different bifurcations.

Note that the main requirement for the scaling law to hold is that the
successive derivatives of $f(x)$ at $x=p$ are zero from order 2 to $k-1$,
independently of the value of $\mu$. In the considered examples, the fact that
$p$ is constant ($p=0$, in particular) is what makes that all the derivatives
from order 2 to $k-1$ are zero (together with the independence of $f^{k}(x)$
on $\mu$). In contrast, when using the stable fixed point $q$, then
$f^{\prime\prime}(q^{*})=0$ but $f^{\prime\prime}(q)\neq 0$ if $q\neq q^{*}$
for $k\geq 3$, and then Eq. (5) and all the subsequent derivations do not
apply. Therefore, the use of a fixed point $p$ that can be stable or unstable
is an important improvement with respect the approach of Ref.
Corral_Alseda_Sardanyes .

Figure 1:  (a) $x_{\ell}$ versus $\mu$ at different values of the time $\ell$
for the discrete dynamical system given by Eq. (10) with $k=3$, yielding the
normal form of the discrete (supercritical) pitchfork bifurcation. Different
initial conditions are taken, uniformly distributed between $q+0.1$ and
$q+0.6$ for $q\geq 0$ and between $q-0.6$ and $q-0.1$ for $q\leq 0$. (b) A
part of the data (for clarity sake) with rescaled axes, representing
$\sqrt{2\ell}\,x_{\ell}$ versus $y=\ell\mu$, as derived from Eq. (8). The
function $\sqrt{G(2y)}$ becomes apparent with the rescaling, and provides the
“microscopic” smooth form of the pitchfork bifurcation (in contrast to the
usual well-known, sharp “macroscopic” picture).

Figure 2:  (a) $x_{\ell}$ versus $\mu$ at different values of the time $\ell$
for the discrete dynamical system given by Eq. (10) with $k=4$. Initial
conditions taken in the same way as in previous figure. (b) A part of the data
(for clarity sake) with rescaled axes, representing
$\sqrt[3]{3\ell}\,x_{\ell}$ versus $y=\ell\mu$, yielding the function
$\sqrt[3]{G(3y)}$. (c) Same data with the more general rescaling obtained from
Eq. (7), displaying ${(k-1)\ell}\,x_{\ell}^{k-1}$ versus
$(k-1)y=(k-1)\ell(\mu-\mu^{*})$, together with other maps, given by $k=2$
(transcritical), $k=3$ (pitchfork, previous figure), and $k=5$, and $f(x)=\mu
x(1-x)$ (logistic, corresponding to $k=2$). Notice that, in some sense, we are
transforming all bifurcations diagrams with $k\geq 3$ into the one of the
transcritical bifurcation (for which $k=2$ and then the exponent $k-1$ does
not introduce any changes to it). $\mu^{*}=1$ for the logistic and $\mu^{*}=0$
for the rest.

##  IV Finite-time scaling law for continuous dynamical systems and
unification with the discrete case

###  IV.1 Theoretical derivation

In this section we consider one-dimensional continuous dynamical systems given
by

$x^{\prime}(t)=g(x),$ (11)

where $x$ is real, $g(x)$ is a real function, and $t$ denotes continuous time.
As in the discrete case, it is assumed that the system has a (stable or
unstable) fixed point at $x=p$, with $g(p)=0$ (and $g^{\prime}(p)<0$ or
$g^{\prime}(p)>0$, respectively, with a bifurcation taking place at
$g^{\prime}(p)=0$) Strogatz_book , and that the initial condition, $x(0)$,
belongs to the basin of attraction of the fixed point. Also, it is implicit
that the function $g(x)$ depends on a bifurcation parameter $\mu$ whose change
may change the value of $p$.

In order to study the transient behavior of the continuous dynamical system,
we may approximate this system by a discrete one, writing Eq. (11) as
$x(t+\Delta)\simeq x(t)+g(x(t))\Delta$ with $\Delta>0$ sufficiently small.
Therefore, we consider the discrete system given by Eq. (1),
$x_{n+1}=f(x_{n})$, with $n$ the discrete time and

$f(x)=x+\Delta g(x).$ (12)

Using a one-step Euler approach numeric , with time step $\Delta$, one can
show that $x_{n}$ is an approximation of the continuous system’s solution, in
the sense of $x_{n}\simeq x(t_{n})$ for $t_{n}=\Delta n$, taking into account
that the initial condition has to fulfill $x_{0}=x(0)$.

Besides, it is straightforward to see that both systems have the same fixed
points, and, due to the fact that $\Delta$ is small enough, the stability (and
the instability) of the discrete system is the same as that of the continuous
one. Indeed, if $g(p)=0$ with $g^{\prime}(p)<0$, then $f(p)=p$ and, for
$\Delta$ small, $|f^{\prime}(p)|=1-\Delta|g^{\prime}(p)|<1$ (but still close
to one, and thus we avoid the case $f^{\prime}(p)\simeq-1$), which are the
conditions for existence and linear stability of fixed points in discrete
dynamical systems Strogatz_book  (on the other side, if $g^{\prime}(p)>0$,
then $f^{\prime}(p)>1$). In addition, $g^{\prime}(p^{*})=0$ implies
$f^{\prime}(p^{*})=1$ and therefore, if the continuous system has a
bifurcation, the discrete system has the same bifurcation, with the same
bifurcation point $\mu=\mu^{*}$. In summary, the whole bifurcation diagram of
the continuous system is inherited by the discrete system. In this way, a
continuous dynamical system can be replaced by a discrete one defined by a
close-to-linear map $f(x)$ (in this sense, continuous dynamical systems can be
considered a special case of discrete ones).

Due to the fact that the continuous and discrete dynamical systems (11) and
(1) are equivalent when the relation (12) holds, the continuous one also
fulfills the scaling law (7)-(8), with the same scaling function $G(y)$ given
by Eq. (4). Naturally, for the continuous system it is more convenient to
express the scaling law in terms of the variables of the continuous system.
Thus, we substitute $f^{\prime}(p)=1+hg^{\prime}(p)$ in
$y=\ell(f^{\prime}(p)-1)$ [Eq. (6)] to get

$y=tM=tg^{\prime}(p),$ (13)

with $t=\ell\Delta$ and $M=g^{\prime}(p)$. Also,
$P^{*}\ell=-tg^{(k)}(p^{*})/k!=C^{*}t$, where we have defined
$C^{*}=-g^{(k)}(p^{*})/k!$ (then, $M=(m-1)/\Delta$ and $C^{*}=P^{*}/\Delta$).
It turns out to be, in the same way as for $f(x)$, that $k$ (integer) is the
lowest order for which $g^{(k)}(p)\neq 0$, with $k\geq 2$.

Therefore, substituting the previous relations into Eq. (8), the scaling law
for the continuous dynamical system (11) can be written as,

$x(t)-p=(\pm 1)^{k}\left({\frac{1}{(k-1)C^{*}t}\,G((k-1)y)}\right)^{1/(k-1)},$
(14)

with $y$ given by Eq. (13) and $G()$ given by Eq. (4). Again, the scaling law
is valid for large $t$, and so for $M=g^{\prime}(p)$ close to zero (as $y=tM$
is finite), and then, if there is a bifurcation, the scaling law turns out to
be valid close to the bifurcation point, given by $M=M^{*}=0$.

###  IV.2 Numerical verification in the continuous case

The validity of the finite-time scaling law for continuous bifurcations, Eq.
(14), is tested in Fig. 3 using the continuous system given by $g(x)=\mu
x-x^{k}$ with $k=5$. Figure 3(a) displays $x(t)$ versus $\mu$ for different
values of $t$ and different initial conditions. Note that, as in the discrete
case, $p=0$, with $M=g^{\prime}(p)=\mu$ (thus, the bifurcation takes place at
$\mu=0$), $y=\mu t$, and $C^{*}=1$. Figure 3(b) shows the results obtained
when the more general rescaling analogous to the one in Eq. (7) (and also
contained in Eq. (14)) is applied, including other bifurcations given by the
$g(x)$ above with different values of $k$. The vertical axis is logarithmic,
to test the validity of the scaling law when $x(t)$ is very close to zero. As
in the discrete case (previous figures), the data collapse is nearly perfect,
which supports the good description provided by the scaling law (14) in the
continuous case at finite times.

Figure 3:  (a) $x(t)$ versus $\mu$ at different times $t$ and using different
initial conditions for the continuous dynamical system given by $g(x)=\mu
x-x^{5}$ (thus, $k=5$). Initial conditions are as in the previous figures. (b)
Some part of the data with rescaled axes, displaying ${(k-1)t}\,x(t)^{k-1}$
versus $(k-1)y=(k-1)tM$, together with other continuous systems with $g(x)=\mu
x-x^{k}$, given by $k=2$ (transcritical), $k=3$ (pitchfork), and $k=4$. Note
that the vertical axis is logarithmic.

###  IV.3 Unification of the discrete and continuous scaling laws

We can easily unify the scaling laws for discrete (1) and for continuous (11)
dynamical systems, Eqs. (8) and (14), respectively. These two equations can be
integrated into a single one, writing it in the form equivalent to Eq. (7),
which leads to

$\left(x(\tau)-p\right)^{k-1}=\frac{1}{(k-1)K^{*}\tau}\,G((k-1)y),$ (15)

where we have introduced a generic time $\tau$ that can be discrete,
$\tau=\ell$, or continuous, $\tau=t$, whereas $x(\tau)$ denotes the solution
at time $\tau$ in both cases. Unifying also the notation for the function
defining the dynamics in Eqs. (1) and (11) as $x(\tau+1)=h(x(\tau))$ and
$dx/d\tau=h(x)$, respectively, we can write

$K^{*}=-\frac{1}{k!}\frac{d^{k}h(p^{*})}{dx^{k}}$

(with $k\geq 2$ the lowest order of the derivative for which $h^{(k)}(p)$ is
not zero) and

$y=\tau(\mathcal{M}-\mathcal{M}^{*}),$ (16)

with

$\mathcal{M}=h^{\prime}(p)$

playing the role of the (alternative) natural bifurcation parameter, with
bifurcation point given by $\mathcal{M}^{*}=h^{\prime}(p^{*})$ (equal to 1 in
the discrete case and to 0 in the continuous case). As always, the scaling
function $G$ is given by Eq. (4).

##  V Scaling law for the distance to the attractor

###  V.1 Finite-time scaling for the approach to the attractor

We have already applied the scaling laws, unified in Eq. (15), to the case in
which the fixed point $p$ can be stable or unstable, corresponding to a change
in stability when crossing a bifurcation, with the only requirement, in
practice, that $p$ is constant (in order that the derivatives of $h(x)$ up to
order $k-1$ are zero at $p$, excluding order one). As at the bifurcation point
the fixed point $p$ changes from stable to unstable (or vice versa), the
resulting $x(\tau)-p$ changes from nearly zero to a different value, as seen
in the previous figures. This holds, for instance, for the bifurcations given
by the $f(x)$ in Eq. (10) and their continuous counterpart (transcritical,
pitchfork…), but not for the saddle-node bifurcation.

However, in addition, the scaling law is also valid when $p$ is replaced by a
stable fixed point $q$ (departing $p$ when this becomes unstable). In such a
case, $x(\tau)-q$ quantifies the “distance” to the attractor and has always to
approach zero; this case is contained also in the scaling law, Eq. (15). But,
despite of the fact that the scaling law is the same, the conditions for its
validity are different. We will have two main subcases. Either $q$ is
constant, as given, e.g., by Eq. (10) for $\mu<0$, where the system will be at
one side of the bifurcation (this constitutes, obviously, a subset of the
previous case given by $p$), or, in the second possibility, $k=2$, and $q$
will represent two different stable fixed points, one at each side of the
bifurcation (in the case in which the bifurcation has stable fixed points at
both sides). The discrete version of this was studied in Ref.
Corral_Alseda_Sardanyes .

In this second subcase the unified scaling law Eq. (15) simplifies to

$x(\tau)-q=\frac{1}{K_{2}^{*}\tau}\,\tilde{G}(z),$ (17)

with $z=\tau(h^{\prime}(q)-h^{\prime}(q^{*}))\leq 0$, from Eq. (16), and
$K_{2}^{*}=-h^{\prime\prime}(q^{*})/2$. The subscript 2 indicates that the
result is only valid for $k=2$, i.e., for the transcritical and the saddle-
node bifurcations. The “new” scaling function is given by

$\tilde{G}(z)=\frac{|z|}{e^{|z|}-1},$

and verifies in fact $\tilde{G}(z)=G(z)$ for $z\leq 0$ ($z$ is defined to be
in that range, where $z=y$), but $\tilde{G}(y)\neq G(y)$ when $y>0$. The
different notation for this scaling function is in order to stress this.

In summary, when the stable fixed point $q$ is constant (at one side of the
bifurcation), the scaling law (15) is perfectly valid, just replacing $p$ by
$q$; and in order to stress that $y$ can only take negative values, we may
replace $G(y)$ by $\tilde{G}(z)$. This holds for the example defined by Eq.
(10), but not for the saddle-node bifurcation. In contrast, when $k=2$
(transcritical and saddle-node bifurcations) we have that the scaling law (15)
transforms into a simplified form, given by Eq. (17). Nevertheless, for the
transcritical bifurcation both subcases are the same.

On the other hand, when $q$ is not constant and $k\geq 3$, the scaling law
(17) does not apply but one can straightforwardly compute the distance to the
attractor if one is able to calculate the value of the attractor (as a
function of the bifurcation parameter), using that the scaling law is not
valid for the attractor $q$ but for another fixed point $p$ if this keeps a
constant value. Obviously, $x(\tau)-q=x(\tau)-p+p-q$, where one just needs to
substitute the scaling law Eq. (15) and the dependence of $q$ with $\mu$ (and
therefore with $\mathcal{M}$). In the particular cases given by Eq. (10), with
$k\geq 2$, and extending them also to continuous dynamical systems, one finds
that $q=0$ for $\mu\leq 0$ and $q=(\pm 1)^{k}\sqrt[k-1]{\mu}$ for $\mu\geq 0$,
leading to

$x(\tau)-q=\frac{(\pm
1)^{k}}{\sqrt[k-1]{(k-1)\tau\,}}\,\left[\sqrt[k-1]{G((k-1)\tau\mu)}-\sqrt[k-1]{(k-1)\tau\mu\,}\,\Theta((k-1)\tau\mu)\right],$
(18)

when $x(0)$ is above the unique attractor $q$ for $k$ even and $|x(0)|$ is
above the absolute value $|q|$ of the two attractors for $k$ odd, with
$\Theta()$ being the Heaviside (unit step) function. Note that the previous
equation constitutes also a scaling law, and it only agrees with Ref.
Corral_Alseda_Sardanyes  for $k=2$, for which $z=-|y|$ and we recover Eq.
(17). This is illustrated in Fig. 4, where the shape of the scaling function
$\tilde{G}(y)$ emerges from the data collapse in the cases in which the
conditions for the validity of the scaling law are fulfilled (in the opposite
case, $k\geq 3$ and $\mu>0$, one can observe how the scaling law is not
fulfilled). Table 2 provides a summary of the rescaling necessary for each
variable in order that the scaling law is fulfilled.

Figure 4:  Distance to the attractor, rescaled, versus distance of the natural bifurcation parameter to the bifurcation point, also rescaled, at different values of the time for the discrete (D) dynamical systems contained in Eq. (10), for the saddle-node bifurcation, given by $f(x)=\mu+x-x^{2}$, and for their continuous (C) counterparts. In concrete, the horizontal axis displays $(k-1)y$, with $y=\tau(\mu-\mu^{*})$ for all cases except for the saddle node, for which $y=-2\tau\sqrt{\mu}$ (we have omited this minus sign for the sake of visualization). $\mu^{*}=0$ in all cases, except for the logistic one. Initial conditions as in previous figures.  Table 2:  Rescaled axes for $\mu$ (horizontal axis) and $(x-p)$ (vertical axis) in order that the scaling laws hold. The two columns on the left refer to particular rescalings for each bifurcation, in which the resulting scaling functions depend on $k$ (tipically as a power of $G()$ with exponent $1/(k-1)$). The two columns on the right lead to a universal scaling for different bifurcations. For the case of Eq. (10) with $k\geq 3$, and its continuous counterparts, the rescaling on the right leads to the scaling function $G()$ when the fixed point is constant, denoted here as $p_{ct}$ (if an attractive fixed point $q$ is used instead, a scaling law still holds, but with different scaling functions, see Eq. (18)). For the saddle-node bifurcation, given by $f(x)=\mu+x-x^{2}$ or $g(x)=\mu-x^{2}$, the scaling law holds for the distance to the attractor $q$, yielding $\tilde{G}$. The transcritical and the logistic (which has also a transcritical bifurcation) fulfill both scaling laws.  bifurcation | horizontal | vertical …. | horizontal | vertical …..
---|---|---|---|---
trans | $\tau\mu$ | $\tau(x-p)$ | $\tau\mu$ | $\tau(x-p)\phantom{{}^{k-1}}$
logistic | $\tau(\mu-1)$ | $\tau(x-p)$ | $\tau(\mu-1)$ | $\tau(x-p)\phantom{{}^{k-1}}$
pitchfork | $\tau\mu$ | $\sqrt{2\tau}\,(x-p_{ct})$ | $2\tau\mu$ | $2\tau(x-p_{ct})^{2}\phantom{{}^{-1}}$
Eq. (10) | $\tau\mu$ | $\sqrt[k-1]{(k-1)\tau}\,(x-p_{ct})$ | $(k-1)\tau\mu$ |  $(k-1)\tau(x-p_{ct})^{k-1}$
saddle-n | $4\tau^{2}\mu$ | $\tau(x-q)$ | $2\tau\sqrt{\mu}$ | $\tau(x-q)\phantom{{}^{k-1}}$
Eq. (19) | $4\tau^{2}\mu$ | $\tau(x-q)$ | $\pm 2\tau\sqrt{|\mu|}$ | $\tau(x-q)\phantom{{}^{k-1}}$

###  V.2 Illustration of validity of the scaling law beyond its limits of
validity

We now analyze the approach to the attractor in a case in which there is not a
strict bifurcation. Let us consider a discrete dynamical system given by

$f(x)=\mbox{sign}(x+\mu)\left(\frac{\sqrt{1+4|x+\mu|}-1}{2}\right),$ (19)

with $\mbox{sign}()$ giving the sign of its argument (being zero when this is
zero). This apparently complicated map arises as the inverse of the very
simple map $f(x)=-\mu+x+\mbox{sign}(x)x^{2}$, where all the fixed points, for
different $\mu$, except $\mu=0$, are linearly unstable (the reason for writing
$\mbox{sign}(x)x^{2}$ and not $x^{3}$ in $f(x)$ is in order to facilitate its
inversion). Therefore, in the case of the map given by Eq. (19), the
“bifurcation diagram” consists of just one linearly stable fixed point (except
for $\mu=0$). This is given by $p=q=-\sqrt{-\mu}$ for $\mu<0$ and
$p=q=\sqrt{\mu}$ for $\mu>0$, (and $p=0$, which is marginal for
$\mu=\mu^{*}=0$). Thus, nothing bifurcates at $\mu=0$ but still we have a
qualitative different behavior just at that single point.

We test the scaling law for such a dynamical system, for which
$y\simeq-2\tau\sqrt{|\mu|}$ (for $\mu$ small) and $K=\mbox{sign}(\mu)$ (except
for $\mu=0$, for which the second derivative at the fixed point is not
defined, and thus there is not a value of $K^{*}$ that we can use). Figure
5(a) displays the solution $x(\tau)$ for different values of $\mu$ at several
times $\tau$, showing how $x(\tau)$ approaches the attractor $q$ for different
values of $\mu$ as $\tau$ increases. Figure 5(b) represents the corresponding
rescaled distance to the attractor as dictated by Eq. (17) (remember, a
particular case of Eq. (15), but in this case without substituting $K$ by
$K^{*}$). We obtain a very good data collapse, meaning that a scaling law
holds, but notice that the scaling function $\tilde{G}(z)$ does not always
describe the rescaled data.

In fact, the scaling function only works for the cases in which the validity
of the scaling law is fulfilled, which corresponds to when $x(\tau)-q$ and $K$
have the same sign, that is, $x(\tau)>q$ when $\mu>0$ and $x(\tau)<q$ when
$\mu<0$. This validity demonstrates that $\mu=0$ is somehow a critical point,
despite one cannot consider that there is a bifurcation there. It is
noteworthy that a scaling law is still fulfilled in the cases in which our
derivation is not valid (corresponding to a range of $\mu$ and $x(\tau)$
opposite to the one given above), and in particular, it is intriguing that the
maximum distance to the attractor $q$ is not attained at the critical point
$\mu=0$, as described by $\tilde{G}(z)$, see Fig. 5(b). From this figure, we
can numerically obtain the value of $\mu$ for which the convergence to the
attractor is slower, which is given by $2\tau\sqrt{|\mu_{\text{slow}}|}\simeq
2.4$, and thus, $\mu_{\text{slow}}\simeq\pm 1.4/\tau^{2}$. In summary, we have
explained one situation in which the scaling law is valid but the scaling
function is not (but still it is valid for a limited range of $\mu$, as
dictated by our analytical calculations).

Figure 5:  (a) $x_{\ell}$ versus $\mu$ at different values of time $\ell$ for
the discrete dynamical system given by Eq. (19). For increasing times, the
solution $x_{\ell}$ approaches $q=\mbox{sign}(\mu)\sqrt{|\mu|}$ (b) Distance
to the attractor with rescaled axes (reducing the number of data points, for
clarity sake). For the sake of visualization, we change the sign of $y$ when
$\mu>0$, i.e., $y=2\tau\sqrt{\mu}$ in this case. A finite-time scaling law is
fulfilled, but the scaling function is only valid (as expected) for a half of
the parameter space. The maximum distance is attained at about 2.4 in the
horizontal axis.

##  VI Conclusions

We have shown how a very important concept in statistical physics, finite-size
scaling, can be used to illuminate the finite-time dynamics in dynamical
systems close to a bifurcation point. Our results extend the previous ones
from Ref. Corral_Alseda_Sardanyes  in two directions: beyond the transcritical
and saddle-node bifurcations and beyond discrete dynamical systems to the
continuous (differential equations) case. Notably, in contrast to Ref.
Corral_Alseda_Sardanyes , our approach does not need to calculate the value of
a stable fixed point, and can proceed with the knowledge of a unstable one
(which in many cases can be simpler to know). The key result is the finite-
time scaling law, Eq. (15), together with the definition of the rescaled
distance $y$ to the bifurcation point, Eq. (16). It is striking the
universality of the scaling function $G(y)$, Eq. (4), which arises (raised to
some power if $k\geq 3$) in all finite-time scaling laws that we obtain.

The parallelisms between bifurcations in dynamical systems and phase
transitions in statistical physics are remarkable, and are summarized in Table
1. In this way, we conclude that the sharp change of behavior that
characterizes bifurcations at their bifurcation point at infinite time arises
in the same way as phase transitions arise in the infinite system-size limit
(known as the thermodynamic limit). The way in which both asymptotic behaviors
arise is through a scaling law, with a smooth scaling function for which the
infinite limit of one of its variables (time or system size) brings the
asymptotic behavior of the scaling function to the transition point
(bifurcation point or critical point). Then, the “critical” exponents
describing the transition are obviously contained in the scaling law and in
the scaling function.

##  VII Acknowledgements

We are grateful to Aina Cornellà, who participated in the early steps of this
work, and to J. Sardanyés, for discussions and a critical reading. The Centre
de Recerca Matemàtica is supported by the CERCA Programme of the Generalitat
de Catalunya as well as by the Spanish State Research Agency (AEI) through the
Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence
in R&D (CEX2020-001084-M). The research of the author is supported by the
projects PGC-FIS2018-099629-B-I00 and PID2021-125979OB-I00, AEI.

##  VIII Appendix: Critical slowing down

###  VIII.1 Power-law decay at the bifurcation point

Naturally, the previous finite-time scaling law (which has two “versions”, one
for the bifurcation profile and another for the distance to the attractor)
include the case in which the system is right at the bifurcation point, $y=0$
(i.e., $\mathcal{M}=\mathcal{M}^{*}$). As $G(0)=1$, we obtain from Eq. (15) a
power-law decay with $\tau$ towards the attractor, given by

$x^{*}(\tau)-q^{*}=\frac{(\pm
1)^{k}}{\left[{(k-1)K^{*}\tau}\right]^{1/(k-1)}},$ (20)

which is in fact the generalization (including both discrete and continuous
systems) of one of the cases in Eq. (9), where we have replaced $p$ by $q^{*}$
as the decay indicates that the fixed point at $\mathcal{M}^{*}$, despite not
being linearly stable, is indeed an attractor. The asterisk in $x^{*}(\tau)$
stresses that we are right at the bifurcation in the parameter space, and the
exponent $1/(k-1)$ quantifies that the higher the value of $k$, the slower the
power-law decay. This equation is valid for the transcritical and pitchfork
bifurcations, as represented by Eq. (10), including higher values of $k$, as
well as for their continuous counterparts, and also for the saddle-node
bifurcation, replacing $p$ by $q$ in Eq. (15). In concrete,
$x^{*}(\tau)-q^{*}\propto 1/\tau$ for the transcritical and saddle-node
bifurcations, $|x^{*}(\tau)-q^{*}|\propto 1/\sqrt{\tau}$ for the pitchfork
case, $x^{*}(\tau)-q^{*}\propto 1/\sqrt[3]{\tau}$ for $k=4$ and so on. A
power-law decay in time when the system is in a critical point is referred to
in physics as critical slowing down Strogatz_book . The fulfillment of the
scaling in the particular case $y=0$ in Figs. 1(b), 2(b)-(c), 3(b), 4, and
5(b) provides visual verification of this critical slowing down.

###  VIII.2 Divergence of the characteristic decay time

The scaling law (15), through the rescaled variable $y$, Eq. (16), relates
values of the time $\tau$, the bifurcation parameter $\mu$ (or $\mathcal{M}$),
and the rescaled distance of the solution to the fixed point that are
“corresponding”, in the sense that they yield the same value of the function
$G()$, and thus, these “corresponding” parameters and values are uniquely
described by the value of $G()$. But the law can also be used to compare the
dynamics at fixed $\mu$ for different $\tau$. Outside the bifurcation point,
we consider first with the case $y<0$. Then, $G(u)=\tilde{G}(u)\rightarrow
ue^{u}$ when
$u=(k-1)y=(k-1)\tau(\mathcal{M}-\mathcal{M}^{*})\rightarrow-\infty$, and
substituting in Eq. (15) we obtain

$x(\tau)-q=(\pm
1)^{k}\sqrt[k-1]{\frac{|\mathcal{M}-\mathcal{M}^{*}|}{K^{*}}}\,e^{-\tau/\tau_{c}}\,\mbox{
for }\mathcal{M}<\mathcal{M}^{*},$ (21)

thus, the decay is exponential, with $\tau_{c}$ a characteristic decay time
given by

$\tau_{c}=\frac{1}{\mathcal{M}^{*}-\mathcal{M}}.$

In the opposite case, $y>0$, we use Eq. (18), which describes the results for
the dynamical system given by Eq. (10). When $u\rightarrow\infty$ we can
substitute

$\sqrt[k-1]{G(u)}-\sqrt[k-1]{u}\rightarrow\frac{1}{k-1}{\sqrt[k-1]{u}\,e^{-u}},$

to arrive to Eq. (21) with $K^{*}=1$ but

$\tau_{c}=\frac{1}{(k-1)(\mathcal{M}-\mathcal{M}^{*})}\mbox{ for
}\mathcal{M}>\mathcal{M}^{*}.$

In this way, the characteristic time $\tau_{c}$ is different at each side of
the bifurcation if $k\geq 3$, and this different behavior is fully contained
in the scaling function. In any case, notice that at the bifurcation point,
$\tau_{c}\rightarrow\infty$ and the exponential decay is no longer valid.
Instead, Eq. (20) applies.

We have shown that $\tau_{c}$ decays hyperbolically in terms of the natural
bifurcation parameter $\mathcal{M}$. For the transcritical bifurcation, the
pitchfork bifurcations, and the cases given by higher values of $k$ arising
from Eq. (10), we have that $\mathcal{M}-\mathcal{M}^{*}=\mu$ (or equal to
$\mu-\mu^{*}$ with $\mu^{*}=1$ for the transcritical bifurcation contained in
the logistic map); so, the decay is the same in terms of $\mu$. However, for
the saddle-node bifurcation
$y=-z=\tau(\mathcal{M}-\mathcal{M}^{*})=2\tau\sqrt{\mu}$, and therefore
$\tau_{c}=1/\sqrt{4\mu}$ in this case. Thus, the critical exponent that
relates $\tau_{c}$ with $\mu$ is $1/2$ in the saddle-node bifurcation, but
turns out to be 1 if one uses the natural bifurcation parameter $\mathcal{M}$.
The map given by Eq. (19) displays the same behavior.

## References

  * (1) S. H. Strogatz.  Nonlinear Dynamics and Chaos.  Perseus Books, Reading, 1994.
  * (2) R. V. Solé. Phase Transitions. Princeton University Press, 2011.
  * (3) T. E. Harris.  The Theory of Branching Processes.  Springer, Berlin, 1963.
  * (4) M. Kimmel and D. E. Axelrod.  Branching Processes in Biology.  Springer-Verlag, New York, 2002.
  * (5) A. Corral and F. Font-Clos.  Criticality and self-organization in branching processes: application to natural hazards.  In M. Aschwanden, editor, Self-Organized Criticality Systems, pages 183–228. Open Academic Press, Berlin, 2013.
  * (6) H. E. Stanley.  Introduction to Phase Transitions and Critical Phenomena.  Oxford University Press, Oxford, 1973.
  * (7) J. M. Yeomans.  Statistical Mechanics of Phase Transitions.  Oxford University Press, New York, 1992.
  * (8) H. E. Stanley.  Scaling, universality, and renormalization: Three pillars of modern critical phenomena.  Rev. Mod. Phys., 71:S358–S366, 1999.
  * (9) R. Garcia-Millan, F. Font-Clos, and A. Corral.  Finite-size scaling of survival probability in branching processes.  Phys. Rev. E, 91:042122, 2015.
  * (10) A. Corral, R. Garcia-Millan, and F. Font-Clos.  Exact derivation of a finite-size scaling law and corrections to scaling in the geometric Galton-Watson process.  PLoS ONE, 11(9):e0161586, 2016.
  * (11) V. Privman.  Finite-size scaling theory.  In V. Privman, editor, Finite Size Scaling and Numerical Simulation of Statistical Systems, pages 1–98. World Scientific, Singapore, 1990\.
  * (12) C. Falcó and A. Corral.  Finite-time scaling for epidemic processes with power-law superspreading events.  arXiv, 2110.10459, 2021.
  * (13) A. Corral, R. Garcia-Millan, N. R. Moloney, and F. Font-Clos.  Phase transition, scaling of moments, and order-parameter distributions in Brownian particles and branching processes with finite-size effects.  Phys. Rev. E, 97:062156, 2018.
  * (14) A. Corral, L. Alsedà, and J. Sardanyès.  Finite-time scaling in local bifurcations.  Sci. Rep., 8(11783), 2018.
  * (15) K. Christensen and N. R. Moloney. Complexity and Criticality. Imperial College Press, London, 2005.
  * (16) A. Corral. Scaling in the timing of extreme events. Chaos. Solit. Fract., 74:99–112, 2015.
  * (17) If one wants to keep track of the sign of $d_{n}$, defining, for instance, $c_{n}=\mbox{sign}(d_{n})/d_{n}^{\alpha}$ for $\alpha$ even, the final results are the same.
  * (18) A. Cornellà. Internal Report, 2021.
  * (19) R. L. Burden and J. D. Faires.  Numerical Analysis.  Brooks Cole, 9 edition, 2010.