11institutetext: Faculdade de Engenharia, Universidade do Porto
11email<EMAIL_ADDRESS>11email:
<EMAIL_ADDRESS>22institutetext: Artificial Intelligence
& Computer Science Lab.(LIACC – member of LASI LA), Universidade do Porto
33institutetext: Fraunhofer AICOS Portugal 44institutetext: INESC TEC,
Portugal

# Kernel Corrector LSTM ††thanks: This work was partially funded by projects
AISym4Med (101095387) supported by Horizon Europe Cluster 1: Health,
ConnectedHealth (n.º 46858), supported by Competitiveness and
Internationalisation Operational Programme (POCI) and Lisbon Regional
Operational Programme (LISBOA 2020), under the PORTUGAL 2020 Partnership
Agreement, through the European Regional Development Fund (ERDF); NextGenAI -
Center for Responsible AI (2022-C05i0102-02), supported by IAPMEI, and also by
FCT plurianual funding for 2020-2023 of LIACC
(UIDB/00027/2020_UIDP/00027/2020) and SONAE IM Labs@FEUP.

Rodrigo Tuna 11 0009-0009-0047-2052 Yassine Baghoussi 1144 0000-0002-1943-1471
Carlos Soares 112233 0000-0003-4549-8917 João Mendes-Moreira 1144
0000-0002-2471-2833

###### Abstract

Forecasting methods are affected by data quality issues in two ways: 1. they
are hard to predict, and 2. they may affect the model negatively when it is
updated with new data. The latter issue is usually addressed by pre-processing
the data to remove those issues. An alternative approach has recently been
proposed, Corrector LSTM (cLSTM), which is a Read & Write Machine Learning
(RW-ML) algorithm that changes the data while learning to improve its
predictions. Despite promising results being reported, cLSTM is
computationally expensive, as it uses a meta-learner to monitor the hidden
states of the LSTM. We propose a new RW-ML algorithm, Kernel Corrector LSTM
(KcLSTM), that replaces the meta-learner of cLSTM with a simpler method:
Kernel Smoothing. We empirically evaluate the forecasting accuracy and the
training time of the new algorithm and compare it with cLSTM and LSTM. Results
indicate that it is able to decrease the training time while maintaining a
competitive forecasting accuracy.

###### Keywords:

Time series forecasting Recurrent Neural Networks Data-Centric AI

##  1 Introduction

In many fields, including energy, healthcare, management, and climate
research, time series forecasting is a crucial task that can be accomplished
using machine learning or statistical methods [8]. As data becomes widely
available, more precise forecasting models are expected. However, data quality
issues like outliers, missing values, and changes in the underlying data
generation process might impact predictive techniques.

Traditional machine learning (ML) models are often considered read-only
models, capable of learning from data but neglecting the feedback loop for
correcting the data during the learning process. This approach, while
efficient in many cases, lacks proper adaptation of preprocessing techniques
and the ML model itself, as the model’s feedback is often overlooked.

To address this limitation, the concept of Read-Write Machine Learning (RW-ML)
has emerged. RW-ML models, such as Corrector LSTM (cLSTM) [1], not only learn
from data but also have the capability to change the data during the learning
process. cLSTM is a time series forecasting method designed to improve
forecasting accuracy by dynamically adjusting the data. It utilizes a meta-
model of the Hidden State Dynamics obtained with SARIMA to detect data quality
issues and employs a greedy heuristic to correct them. cLSTM has demonstrated
superior predictive performance compared to traditional LSTM models. However,
the computational cost associated with the meta-learning component of cLSTM is
significant.

In this paper, we propose a computationally less expensive variant of cLSTM,
named Kernel Corrector LSTM (KcLSTM), which replaces the meta-learner with a
simpler method: Kernel Smoothing. We empirically compare KcLSTM with both
cLSTM and LSTM models. Results reveal that KcLSTM achieves better predictive
performance than LSTM and cLSTM, while also being faster than cLSTM, although
the computational efficiency improvement is not as substantial as expected.

The main contributions of this paper are:

  * •

Introducing a variant of cLSTM, KcLSTM, which is computationally less
expensive while maintaining high predictive accuracy.

  * •

An empirical study comparing KcLSTM with LSTM and cLSTM in terms of predictive
performance and training time.

This paper is structured as follows: we first provide an overview of the
state-of-the-art forecasting method, LSTM. Then, we delve into the concept of
RW-ML and its significance in time series forecasting. Next, we introduce the
proposed algorithm, KcLSTM. Finally, we describe the experimental setup,
present the results, and discuss their implications.

##  2 Related Work

In this section, we first present the Long Short-Term Memory. The algorithm
that our proposed algorithm is built on and the one it will be compared to. We
then define Data-Centric AI and provide examples of Data-Centric models built
for time series forecasting.

###  2.1 LSTM

Figure 1: Cell unit of the LSTM recurent neural network [6].

The Long Short-Term Memory (LSTM) [16] is a Recurrent Neural Network (RNN),
that can capture long-term dependencies in the input and it is used for
processing sequential data. RNNs differ from feed-forward networks through
recurrent connections, allowing them to learn from sequential data. Back
Propagation is applied to RNNs by taking advantage of the fact that for every
recurrent network, there exists an equivalent feed-forward network with
identical behavior for a finite number of steps [21], training it using Back
Propagation Through Time (BPTT) [23].

RNNs have some well-known limitations. First, they have problems capturing
long-term dependencies, being limited to only bridge between 5-10 steps [23].
This occurs because RNNs are sensible to the exploding/vanishing gradient
problem [15]. The LSTM solves this problem through the use of a gating
mechanism.

An LSTM network consists of blocks, with each block containing an input gate,
forget gate, output gate, and memory cell (eqs. 1b, 1c and 1d). The input gate
controls which inputs are relevant; the forget gate learns which information
should be kept in memory; and the output gate controls which information
should be passed to the next block. The information is retained through the
use of two states called the cell state, (eq. 1e), and the hidden state, (eq.
1f). The forward pass concatenates the input with the hidden state from the
last block, while the backward pass derives the error and updates the gates
using the chain rule of derivatives. The gate derivatives are multiplied by
the hidden output to obtain the gradient deltas that update the gates.

$\displaystyle i_{t}$ $\displaystyle=\sigma(W_{i}\cdot h_{t-1}+V_{i}\cdot
x_{t}+b_{i})$ (1a) $\displaystyle o_{t}$ $\displaystyle=\sigma(W_{o}\cdot
h_{t-1}+V_{o}\cdot x_{t}+b_{o})$ (1b) $\displaystyle f_{t}$
$\displaystyle=\sigma(W_{f}\cdot h_{t-1}+V_{f}\cdot x_{t}+b_{f})$ (1c)
$\displaystyle\hat{C_{t}}$ $\displaystyle=tanh(W_{c}\cdot h_{t-1}+V_{c}\cdot
x_{t}+b_{c})$ (1d) $\displaystyle C_{t}$
$\displaystyle=i_{t}\cdot\hat{C_{t}}+f_{t}\cdot C_{t-1}$ (1e) $\displaystyle
h_{t}$ $\displaystyle=o_{t}\cdot tanh(C_{t})$ (1f) $\displaystyle z_{t}$
$\displaystyle=h_{t}$ (1g)

###  2.2 Data-Centric Time Series Forecasting

Anomalies, including outliers, missing values, and changes in the underlying
data generation process can impact predictive tasks. This affects the
predictions of such methods, hindering their performance [17], conversely to
traditional machine learning methods, that build models using a fixed dataset.
In Data-Centric AI [25], the focus is on the data. Data quality is increased
to improve the performance of AI models.

The exploration into machine learning models capable of learning and
correcting data has been a topic of interest in various studies. Both [22] and
[20] delve into this concept, with [22] focusing on the potential of ML models
to memorize sensitive information while [20] emphasize the importance of model
interpretability and safety. Additionally, authors in [4] further underscore
the significance of data quality in enhancing model performance, advocating
for a data-centric approach. Providing a broader perspective, the work in [14]
discusses the role of probabilistic modeling in understanding learning and
uncertainty in machine learning.

Moving to neural network models, authors in [11] discuss highly interconnected
networks for associative memory and optimization, with a focus on learning and
adaptation. Moreover, [10] propose a model for neural networks that learn
temporal sequences through selection, employing synaptic triads and a local
Hebbian learning rule. Furthermore, [9] introduce predictive-corrective
networks for action detection in videos, which utilize top-down predictions
and bottom-up observations for adaptive computation and simplified learning.
These models collectively demonstrate the potential of neural networks to
learn and correct data across various applications.

Similarly, recurrent neural network (RNN) models have been developed to
address the challenge of learning and correcting data. In [13], an attempt is
made to introduce a learning algorithm for the recurrent random network model,
employing gradient descent of a quadratic error function. Later, authors in
[2] propose a recurrent network architecture for modeling dynamical systems,
which can learn from multiple temporal patterns and cope with sparse data.
More recently, the research in [5] demonstrates that tree-structured recursive
neural networks can learn logical semantics, including entailment and
contradiction.

In the context of time series forecasting, some data-centric approaches have
been employed. In dLSTM [19], the authors train the model on non-anomalous
data and use the predictive errors to detect anomalies. The deviation from the
normal state is measured through delayed prediction errors. The normal state
can then be restored from several candidate values. Following the idea of
using the prediction errors to improve the quality of the data Pastprop was
introduced [3]. The responsibility for the training error is shared between
the model parameters and the training data. The backpropagation of the
derivatives is applied to the input, indicating the part of the input that
caused the training error.

###  2.3 Corrector LSTM

cLSTM [1] is an architecture that improves its predictive performance by
reconstructing the data of the model. The architecture of the algorithm is
based on the LSTM and a data correction component. This data correction
component uses a meta-learner, SARIMA, to identify problems in the hidden
states of an LSTM model. This is achieved by predicting the hidden states
using SARIMA and if the difference between the predicted and the real hidden
states is over a certain threshold they are considered anomalous. The
anomalies detected in the hidden states are assumed to be caused by the data
which is then reconstructed. The reconstruction of the data points is such
that the difference between the predicted and the real hidden states falls
under a certain threshold. The authors showed that analyzing the Hidden State
Dynamics [24] of an LSTM can be used to detect anomalies in the training data
and consequently improve the forecasting performance of the model. However,
the data correction relies on a meta-learner which makes the algorithm
computationally expensive.

##  3 Kernel Corrector LSTM

The architecture of the Kernel Corrector LSTM (KcLSTM) is the same as the
cLSTM architecture, and the meta-learner used to detect problems in the
learning is substituted by a simpler approach, kernel smoothing.

###  3.1 Training

The KcLSTM utilizes the hidden states learned during the training process to
find and correct data points of the series. The training of the KcLSTM is
divided into three distinct phases. The first phase consists of training the
data on a standard LSTM. This allows the hidden states to capture the
information of the time series and to be indicative of problems in the data.
We then perform the correction, which is comprised of a detection and a
correction component. These two components find and correct errors in the data
respectively, this phase is thoroughly explained in section 3.2. Finally, the
LSTM is trained on the new data, learning a corrected time series, that can
improve the predictions of the model.

###  3.2 Data Correction

The Data Correction phase of the algorithm is divided into two different
components: the correction and the detection. These two components aim to find
data points that worsen the learning of the model and change the data so that
the learning process is improved and a better model is obtained. Each phase
has a threshold $\delta_{d}$ and $\delta_{c}$. The hidden states of the last
iteration of the first training phase, $H={h_{0},...,h_{n}}$ are used to find
errors in the training data. cLSTM uses a meta-model that is computationally
expensive to compute; our goal is to assess if a simpler method can obtain
competitive results with less cost; the method selected for this purpose is
Kernel Smoothing because states are estimated rather than predicted which
makes it computationally. A new set of estimated hidden states
$H^{\prime}={h^{\prime}_{0},...,h^{\prime}_{n}}$ is calculated using Gaussian
Kernel Smoothing of $H$ as described in eq. 2.

$h^{\prime}_{i}=\frac{\sum_{j\in[i-W/2,i+W/2],i\neq
j}h_{j}*K(h_{i},h_{j})}{\sum_{j\in[i-W/2,i+W/2],i\neq j}K(h_{i},h_{j})}$ (2)

Where $K(h_{i},h_{j})$ is:

$K(h_{i},h_{j})=e^{\frac{\lVert h_{i}-h_{j}\rVert^{2}}{2\sigma^{2}}}$ (3)

Figure 2: Gaussian Smoothing of the Hidden States, represented as a series.

The goal of error detection is to discriminate between data points that need
reconstruction and those that do not. A point needs reconstruction if the
Dynamic Time Warp similarity between the hidden state from which the point
originated $h_{i}$ and the corresponding estimated hidden state
$h^{\prime}_{i}$ is greater than a given threshold $\delta_{d}$. This relation
is depicted in eq. 4.

$DTW(h^{\prime}_{i},h_{i})>\delta_{d}$ (4)

In the error correction, we reconstruct the detected points such that the
Dynamic Time Warp similarity of the hidden state and estimated hidden state is
less or equal to a given threshold $\delta_{c}$ eq. 5. Early stopping is
employed, and if a maximum number of iterations is reached, the original value
for the point is restored.

$DTW(h^{\prime}_{i},h_{i})\leq\delta_{c}$ (5)

##  4 Experimental Setup

The goals of the empirical validation are to investigate if the proposed
algorithm is faster than the original one, without a significant decrease in
forecasting accuracy.

A straightforward holdout method was used to estimate forecasting performance,
used when there is a temporal dependency in the dataset [7]. The model is
trained on the first $s$ samples and assessed on the succeeding $n-s$ samples.
The data used for evaluation is always the original one. Using corrected data
for the evaluation would likely lead to inadequate optimistic estimates of the
forecasting performance of the corresponding method.

The hyperparameters of LSTM and KcLSTM are chosen using hyperparameter tuning
using grid search. The learning rate was varied between: 0.0001, 0.001, 0.01,
0.1; and the batch size was varied between: 1, 2, 4, 8. For cLSTM we do not
perform hyper-parameter tuning due to the high computational cost. Instead, we
use the results described in the original paper [1]. The thresholds for KcLSTM
are fixed with values of 0.6 for the detection and 0.5 for the correction. The
thresholds for cLSTM are described in the original paper [1], 0.6 for the
detection and 0.2 for the correction. We chose to maintain the same detection
threshold and increased the correction threshold. The kernel estimates of the
hidden states are smoother; thus, a small correction threshold would
significantly alter the hidden states, and the information learned in the
previous phase would be lost.

###  4.1 Datasets

Table 1: Statistical description of the dataset. | Monthly
---|---
Timeseries | 200
Average Length | 366
Mean | 4222
Standard Deviation | 1160

We have used the M4 Competition Dataset [18] comprising six subsets. From one
subset, Monthly, we evaluate the performance of the algorithms on the first
199 time series. To evaluate the time taken to train the models, we use the
first 20 time series of that subset.

###  4.2 Evaluation Metrics

This study focuses on both the predictive performance of the algorithm as well
as its training time. To quantify the error of forecasts, we focus on the Mean
Absolute Scaled Error (MASE) in Eq 6 because it allows for the averaging of
results across different time series as opposed to the Rooted Mean Squared
Error (RMSE). The MASE measures the appropriateness of a forecast against the
naive forecast of predicting the previous value. To assess if the differences
are statistically different, we use the Mariano-Diebold Test [12].

$MASE=\frac{\frac{1}{n-s}\sum_{i=s+1}^{n}{\hat{y}_{i}-y_{i}}}{\frac{1}{n-1}\sum_{i=2}^{n}|y_{i}-y_{i-1}|}$
(6)

The execution time is measured in seconds and the experiments were run on an
Intel(R) Core(TM) i7-1065G7 CPU @ 1.30GHz processor.

##  5 Results

investigate if the proposed algorithm is faster than the original one, without
significantly decreasing forecasting accuracy. To illustrate the usefulness of
the proposed algorithm, we first evaluate its forecasting capabilities,
comparing it with LSTM section 5.2 and cLSTM section 5.1.

###  5.1 Comparison with cLSTM

Results of MASE for the algorithms presented in table 2 indicate KcLSTM
outperforms cLSTM, but this may be explained by the hyper-parameter tuning
that was performed for the KcLSTM and not the cLSTM. As such comparison
between these two algorithms, can not be performed directly.

The Mariano-Diebold Test for cLSTM and KcLSTM resulted in 40 wins for cLSTM,
111 wins for KcLSTM, and 48 draws. This shows an improvement in forecasting
accuracy by substituting the meta-learner with Kernel Smoothing. Again the
uneven conditions do not allow us to reach clear conclusions about these two
methods.

Results for the training time presented in table 2 indicate that KcLSTM is
indeed faster than cLSTM, significantly. Nonetheless, the gain is not as large
as would be expected. Estimating the states with Kernel Smoothing is less
computationally expensive than predicting the states with SARIMA. However,
this is a cruder method that results in estimated states that are farther away
than from the original states when compared with. Consequently, more points
are considered anomalies that will be corrected, which will cause the training
time to increase.

Table 2: Comparison of each algorithm. | Mean | Median | Standard Deviation | Average time (s)
---|---|---|---|---
LSTM | 3.48 | 0.74 | 6.07 | 20.47
cLSTM | 8.77 | 1.04 | 36.96 | 56.15
KcLSTM | 4.64 | 0.83 | 11.96 | 48.77

However, when performing the Mariano-Diebold test to compare LSTM and KcLSTM
at the significance level of $\alpha=0.05$ we get 102 wins for KcLSTM, 50 wins
for LSTM, and 47 draws. We can conclude that KcLSTM is superior to LSTM as it
wins more often, although when it loses it is by a greater margin. This is
confirmed by the values of the standard deviation of MASE for LSTM and KcLSTM
and explains the (apparent) superiority of LSTM when analyzing only the MASE.
We see examples of a series with clear outliers that KcLSTM is able to correct
and as such increase their predictions in fig. 3. Conversely, an example of a
series without outliers made KcLSTM wrongfully alter the data which results in
disastrous predictions in fig. 4. These two examples reflect the different
behaviors mentioned before.

Figure 3: Example where data reconstruction was successful. Figure 4: Example
where the data reconstruction destroyed the data which caused the model not to
capture the information of the series.

###  5.2 Comparison with LSTM

Results of MASE for the algorithms presented in table 2 indicate that LSTM has
an overall better performance than KcLSTM with lower values for both the
median and the mean for the MASE. However, when performing the Mariano-Diebold
test to compare LSTM and KcLSTM at the significance level of $\alpha=0.05$ we
get 102 wins for KcLSTM, 50 wins for LSTM, and 47 draws. We can conclude that
KcLSTM is superior to LSTM as it wins more often, although when it loses it is
by a greater margin. This is confirmed by the values of the standard deviation
of MASE for LSTM and KcLSTM and explains the (apparent) superiority of LSTM
when analyzing only the MASE. We see examples of a series with clear outliers
that KcLSTM is able to correct and as such increase their predictions in fig.
3. Conversely, an example of a series without outliers made KcLSTM wrongfully
alter the data which results in worse predictions in fig. 4. These two
examples reflect the different behaviors mentioned before.

Results for the training time presented in table 2 indicate that KcLSTM is
significantly slower than LSTM. This to be expected as KcLSTM has a data
correction component that is responsible for most of the execution time of the
algorithm.

##  6 Conclusion

The goal of our work is to create a forecasting algorithm that reconstructs
data that is faster than current solutions in the literature. We present a new
algorithm: Kernel Corrector LSTM (KcLSTM). This algorithm alters the training
data to improve its forecasting accuracy. Like in cLSTM, this is done by
analyzing the Hidden States Dynamics and finding anomalies in hidden states to
detect anomalies in data points and consequently correct them. However, the
meta-learner of cLSTM was replaced by the Gaussian Kernel Smoothing of the
hidden states to decrease the training time of cLSTM.

We empirically compare our algorithm with LSTM and cLSTM both in terms of
predictive performance and training time. Results show that KcLSTM obtains a
competitive forecasting accuracy surpassing both the LSTM and cLSTM in the
number of statistically significant wins. However, KcLSTM is more sensitive to
its training data and more prone to making worse forecasts than the baseline,
which caused the average MASE of LSTM to be inferior to the average MASE of
KcLSTM. The measured training times also show that KcLSTM indeed improves on
cLSTM in terms of computational cost, but the margin is smaller than expected
because KcLSTM detects more points as anomalies than cLSTM. The estimated
hidden states by KcLSTM are more distant from the real hidden states than the
predicted states of cLSTM. The empirical study showed that KcLSTM is a faster
algorithm that corrects its training data than cLSTM and that those
corrections improve the forecasts by being superior to LSTM. Future work
comprises the possibility of implementing the algorithm with different
estimators.

## References

  * [1] Baghoussi, Y., Soares, C., Mendes-Moreira, J.: Corrector lstm: Built-in training data correction for improved time series forecasting. In: Proceedings of the 8th SIGKDD International Workshop on Mining and Learning from Time Series–Deep Forecasting: Models, Interpretability, and Applications. pp. 1–8. ACM, Washington DC, USA (August 2022)
  * [2] Bailer-Jones, C., MacKay, D.J.C., Withers, P.J.A.: A recurrent neural network for modelling dynamical systems. Network 9 4, 531–47 (1998), https://api.semanticscholar.org/CorpusID:653765
  * [3] Baptista, A., Baghoussi, Y., Soares, C., Mendes-Moreira, J., Arantes, M.: Pastprop-rnn: improved predictions of the future by correcting the past (2021)
  * [4] Bhowmik, P., Partha, A.S.: A data-centric approach to improve machine learning model’s performance in production. International Journal of Engineering and Advanced Technology (2021), https://api.semanticscholar.org/CorpusID:240328155
  * [5] Bowman, S.R., Potts, C., Manning, C.D.: Recursive neural networks can learn logical semantics. In: Workshop on Continuous Vector Space Models and their Compositionality (2014), https://api.semanticscholar.org/CorpusID:15618372
  * [6] Castro, J., Achanccaray Diaz, P., Sanches, I., Cue La Rosa, L., Nigri Happ, P., Feitosa, R.: Evaluation of recurrent neural networks for crop recognition from multitemporal remote sensing images (11 2017)
  * [7] Cerqueira, V., Torgo, L., Mozetič, I.: Evaluating time series forecasting models: an empirical study on performance estimation methods. Machine Learning 109(11), 1997–2028 (oct 2020). https://doi.org/10.1007/s10994-020-05910-7, https://doi.org/10.1007\%2Fs10994-020-05910-7
  * [8] Cerqueira, V., Torgo, L., Soares, C.: Machine learning vs statistical methods for time series forecasting: Size matters (2019)
  * [9] Dave, A., Russakovsky, O., Ramanan, D.: Predictive-corrective networks for action detection. 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) pp. 2067–2076 (2017), https://api.semanticscholar.org/CorpusID:2466592
  * [10] Dehaene, S., Changeux, J.P., Nadal, J.P.: Neural networks that learn temporal sequences by selection. Proceedings of the National Academy of Sciences of the United States of America 84 9, 2727–31 (1987), https://api.semanticscholar.org/CorpusID:7423734
  * [11] Denker, J.S.: Neural network models of learning and adaptation. Physica D: Nonlinear Phenomena 2, 216–232 (1986), https://api.semanticscholar.org/CorpusID:119988262
  * [12] Diebold, F., Mariano, R.: Comparing predictive accuracy. Journal of Business & Economic Statistics 13(3), 253–63 (1995), https://EconPapers.repec.org/RePEc:bes:jnlbes:v:13:y:1995:i:3:p:253-63
  * [13] Gelenbe, E.: Learning in the recurrent random neural network. Neural Computation 5, 154–164 (1992), https://api.semanticscholar.org/CorpusID:38667978
  * [14] Ghahramani, Z.: Probabilistic machine learning and artificial intelligence. Nature 521, 452–459 (2015), https://api.semanticscholar.org/CorpusID:216356
  * [15] Hochreiter, S.: The vanishing gradient problem during learning recurrent neural nets and problem solutions. International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems 6, 107–116 (04 1998). https://doi.org/10.1142/S0218488598000094
  * [16] Hochreiter, S., Schmidhuber, J.: Long short-term memory. Neural computation 9, 1735–80 (12 1997). https://doi.org/10.1162/neco.1997.9.8.1735
  * [17] Kanarachos, S., Christopoulos, S.R.G., Chroneos, A., Fitzpatrick, M.E.: Detecting anomalies in time series data via a deep learning algorithm combining wavelets, neural networks and hilbert transform. Expert Systems with Applications 85, 292–304 (2017). https://doi.org/https://doi.org/10.1016/j.eswa.2017.04.028, https://www.sciencedirect.com/science/article/pii/S0957417417302737
  * [18] Makridakis, S., Spiliotis, E., Assimakopoulos, V.: The m4 competition: Results, findings, conclusion and way forward. International Journal of Forecasting 34(4), 802–808 (2018). https://doi.org/https://doi.org/10.1016/j.ijforecast.2018.06.001, https://www.sciencedirect.com/science/article/pii/S0169207018300785
  * [19] Maya, S., Ueno, K., Nishikawa, T.: dlstm: a new approach for anomaly detection using deep learning with delayed prediction. International Journal of Data Science and Analytics 8 (09 2019). https://doi.org/10.1007/s41060-019-00186-0
  * [20] Otte, C.: Safe and interpretable machine learning: A methodological review (2013), https://api.semanticscholar.org/CorpusID:56899177
  * [21] Sherstinsky, A.: Fundamentals of recurrent neural network (rnn) and long short-term memory (lstm) network. Physica D: Nonlinear Phenomena 404, 132306 (2020). https://doi.org/https://doi.org/10.1016/j.physd.2019.132306, https://www.sciencedirect.com/science/article/pii/S0167278919305974
  * [22] Song, C., Ristenpart, T., Shmatikov, V.: Machine learning models that remember too much. Proceedings of the 2017 ACM SIGSAC Conference on Computer and Communications Security (2017), https://api.semanticscholar.org/CorpusID:2904063
  * [23] Staudemeyer, R.C., Morris, E.R.: Understanding lstm – a tutorial into long short-term memory recurrent neural networks (2019)
  * [24] Strobelt, H., Gehrmann, S., Pfister, H., Rush, A.M.: Lstmvis: A tool for visual analysis of hidden state dynamics in recurrent neural networks (2017)
  * [25] Zha, D., Bhat, Z.P., Lai, K.H., Yang, F., Jiang, Z., Zhong, S., Hu, X.: Data-centric artificial intelligence: A survey (2023)


[1]Kimberly Ayers

These authors contributed equally to this work.

[1]Department of Mathematics, California State University San Marcos, 333 S.
Twin Oaks Valley Rd., San Marcos, 92096, CA, USA

2]Department of Mathematics and Statistics, Pomona College, 610 College Ave.,
Claremont, 91711, CA, USA

# Noisy fixed points: stability of the invariant distribution of the random
logistic map

<EMAIL_ADDRESS>Ami Radunskaya<EMAIL_ADDRESS>*  [

###### Abstract

The full family of discrete logistic maps has been widely studied both as a
canonical example of the period-doubling route to chaos, and as a model of
natural processes. In this paper we present a study of the stochastic process
described by iterations of the logistic map with a random parameter. In
addition to bringing together previously known results, we present a proof
that the unique invariant measure of the process in the chaotic regime is
asymptotically stable.

###### keywords:

stochastic quadratic map; discrete chaos; Markov process; asymptotic stability

##  1 Introduction

The discrete logistic map is well known as a model of self-limiting population
growth, whose dynamic properties are described in the early review paper by
Robert May [1]. In this highly influential paper (over 6,600 citations as of
November, 2023), May states that, while the fine structure of the bifurcation
diagrams in unimodal maps are fascinating, what is more important in practice
is the statistical properties of these maps. In his opinion, “What is called
for is an effectively stochastic description of the dynamics”, and that these
descriptions of evolving continuous densities will still capture some salient
features of the maps themselves. One way to see beyond the micro-structure
created by evolving parameter values in the deterministic setting is to make
these parameters random variables. To reflect what might happen in practice,
where parameter values representing growth rates might fluctuate over some
specific range, the support of these random variables might be a small
interval. Since May’s remarks appeared, many models of this sort have appeared
in the ecology, biology and economics literature. In [2], the authors address
the Neubert conjecture by looking at a stochastically-perturbed Beverton-Holt
model of population growth. A continuous version of the stochastic logistic
model is studied in [3], where they show that these systems can be used to
describe experimental results on microbial growth. Three different economic
models described by randomly perturbed deterministic maps are studied in [4],
[5] and [6]. These examples are by no means exhaustive: they do give a sense
of the ubiquity and impact of these models. In most cases, the goal is to
first show that the map has a unique invariant measure, since the uniqueness
of this measure is essential if one is to say something about the long-term
behavior of the system.

Along these same lines, Kolmogorov suggested starting with a deterministic
dynamical system, and adding a stochastic perturbation in such a way that the
resulting process has a unique invariant measure. This stochastic version of
the original deterministic system is sometimes called the Kolmogorov measure
of the dynamical system [7].

In this paper we present a study of the stochastic process described by
iterations of the logistic map with a randomly-varying parameter. In addition
to bringing together previously known results, we present a proof that the
unique invariant measure of the process is asymptotically stable.

##  2 Background and Definitions

Given a state space $S$, one can consider either a deterministic or a
stochastic (or random) discrete dynamical system on $S$. A deterministic
dynamical system is given by a function $f:S\rightarrow S$, where, given an
initial state $x_{0}\in S$, one can compute future states recursively:

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

Because, given $x_{n}$, the value of $x_{n+1}$ can be computed exactly, this
is a deterministic process, and the set of iterates,
$x_{0},x_{1},x_{2},\dots=\\{x_{n}\\}_{n=0}^{\infty}$, is called the trajectory
of the initial value $x_{0}$ under the dynamical system given by the function
$f$. A natural question is: given an initial value, $x_{0}$, what is the long-
term behavior of the sequence $x_{0},x_{1},x_{2},\dots$? If the sequence
converges to a finite value, it must be a fixed point, i.e. a value $x^{*}$
that satisfies $x^{*}=f(x^{*})$. So one of the first questions a dynamicist
asks when encountering a dynamical system is: where are the fixed points?
Other questions of interest include determining the stability of fixed points
- whether nearby trajectories converge toward, or away from, fixed points -
and looking at the limiting behavior of non-fixed points.

In a stochastic dynamical system, we consider a sequence of random variables:
$\\{X_{n}\\}_{n=0}^{\infty}$. Given a value for the random variable $X_{n}$ in
the state space, the value of $X_{n+1}$ can not be determined exactly because
the law of evolution is subject to some noise. We are then interested in the
sequence of distributions of these random variables. In particular, we are
interested in the existence of an _invariant distribution_ : a distribution of
$X_{n}$ such that $X_{n+1}$ has the same distribution. If such a distribution
exists, we are also interested in its stability; that is, whether other
distributions converge to this invariant distribution in some manner. We
define the random logistic map and the notions of invariant measures and
stability more precisely below.

In this paper, we consider the state space $S=[0,1]$, and the parameter space
$\Lambda=(3.87,4)$. Let $\mathcal{B}(S)$ denote the Borel sigma-algebra on
$S$. We consider the logistic Markov process

$X_{n+1}=\lambda_{n+1}X_{n}(1-X_{n}),$ (1)

where $\\{\lambda_{i}\\}_{i\geq 1}$ are independent, $\Lambda$-valued random
variables, independent of $X_{n}$ with identical distribution $Q$. The initial
state $X_{0}$ is an $S$-valued random variable, independent of the
$\lambda$’s. Let $A\subset S$. Then $\\{X_{n}\\}$ is an $S$-valued Markov
process with a continuous state space, $S$, and transition probabilities

$\hbox{Prob}\\{X_{1}\in A|X_{0}=x\\}=p(x,A)=Q\left(\frac{A}{x(1-x)}\right).$

In words: the probability that $X_{1}$ is in a given set $A$, given $X_{0}=x$,
is the probability of the set of $\lambda$-values that map $x$ to a point in
$A$. The measure $p$ is called a (one-step) transition probability on the
state space, $S$, and we can write $p(x,A)=\int_{S}{\bf 1}_{A}\,p(x,dy)$.

Similarly, we write the $n$th order transition probability as

$\hbox{Prob}\\{X_{n}\in A|X_{0}=x\\}=p^{n}(x,A).$

Suppose we have a distribution, $\mu$, on the state space, $S$. When we apply
the process (1) to subsets of $S$, we induce a new distribution, $P\mu$ on the
state space as follows. If $A\subset S$, then we define

$P\mu(A)=\int_{S}p(x,A)\mu(dx).$ (2)

So the operator $P$ maps probability measures on $S$ to probability measures
on $S$. We can now define an invariant measure under the operator, $P$, which
corresponds to a fixed point in the deterministic case.

######  Definition 1.

If $P\mu(A)=\mu(A)$ for all $A\in{\cal B}(S)$ then $\mu$ is an invariant
measure under the operator $P$.

We can iterate the operator $P$ a positive number $m$ times, which we denote
by $P^{m}$: if $X_{0}$ has distribution $\mu$ then $X_{m}$ has distribution
$P^{m}\mu$. This leads us the analog to stability in the stochastic context.

######  Definition 2.

A sequence of probability measures $P^{n}(n\geq 1)$ on $S$ converges weakly to
a probability measure $P$ on $S$ if, for any continuous, bounded function
$h:S\mapsto S$,

$\lim_{n\rightarrow\infty}\int_{S}h\,dP^{n}=\int_{S}h\,dP$

In the deterministic case, we know that, if the iterates of $x_{0}$ under a
continuous function, $f$, converge to a point, $x$, then $x$ must be a fixed
point of $f$. This result follows from the continuity of $f$. We can get a
similar result in the stochastic setting, but we first need to define a notion
of continuity. Suppose we have a Markov process on $S$ with transition
probabilities given by $p(x,dy)$.

######  Definition 3.

The transition probability $p(x,dy)$ is weakly continuous or has the Feller
property if, for every continuous and bounded function $h:S\mapsto S$, the
function

$x\mapsto\int_{S}h(y)p(x,dy)$

is continuous.

Using weak continuity we get the following result ([8], Theorem 11.1).

######  Proposition 1.

Suppose $p(x,dy)$ is a transition probability on $S$ with the Feller property.

  1. 1.

If there is an $x\in S$, a sequence of integers $1\leq n_{1}<n_{2}<\dots$ and
a probability measure $\pi_{x}$ on $S$ such that

$\frac{1}{n_{k}}\sum_{m=0}^{n_{k}-1}p^{(m)}(x,dy)\hbox{ converges weakly to
}\pi_{x}\hbox{ as }k\rightarrow\infty$

then $\pi_{x}$ is an invariant measure for $p$.

  2. 2.

If the limit above exists for all $x\in S$ and is independent of $x$ then
$\pi=\pi_{x}$ is the unique invariant probability for $p$.

We can state this in terms of the operator $P$ given in Equation (2). Suppose
a process has the Feller property and we have a measure $\mu$ on $S$. If we
have weak convergence in terms of the average to a measure, $\pi$:

$\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{m=0}^{n-1}P^{m}\mu(A){\rightarrow}\pi(A)$

then $\pi$ is an invariant measure under the operator $P$, i.e.
$P\pi(A)=\pi(A)$ for $A\in{\cal B}(S)$. Furthermore, if Condition 2 holds in
Proposition 1, then this invariant measure is unique.

What we’d like is the asymptotic stability of the (unique) invariant measure,
$\pi$. In other words, for any probability measure $\mu$ on $S$, when is it
the case that

$\sup_{A\in{\cal B}(S)}|P^{n}\mu(A)-\pi(A)|\rightarrow 0\quad\hbox{ as
}n\rightarrow\infty?$

To answer this question, we need the notions of irreducibility and
aperiodicity.

The following definitions of irreducibility and recurrence are from [8], Def.
9.2

######  Definition 4.

Let $\mathcal{B}(S)$ denote the Borel subsets of $S$.

A Markov process $\\{X_{n}\\}_{n=0}^{\infty}$ is Harris irreducible with
reference measure $\nu$ (or $\nu$-irreducible for short) if, for every $x\in
S$ and every $B\in\mathcal{B}(S)$ with $\nu(B)>0$, there is $n\geq 1$ such
that $P^{n}(x,B)>0$.

In words: a Markov process is Harris irreducible if every element in the state
space eventually gets mapped to any set of positive probability.

We will also use a “local” version of irreducibility:

######  Definition 5.

Let $J\subset S$ be a set of positive measure. We say that $\\{X_{n}\\}$ is
$\nu$-irreducible on $J$ if, for all $x\in S$ and $B\in\mathcal{B}(J)$ with
$\nu(B)>0$, then there is an $n\geq 1$ such that $P^{n}(x,B)>0$.

The Krylov-Bogolubov theorem [9] tells us that a stochastic dynamical system
obtained from a continuous function on a compact subset of ${\mathbb{R}}^{d}$
with independent parametric noise always has an invariant distribution:

######  Theorem 1 (Krylov-Bogolubov).

A Feller continuous, stationary Markov process on a compact subset of
${\mathbb{R}}^{d}$ has an invariant distribution.

The theorem applies to more general settings, but this more specific version
is all we need here. Harris irreducibility can then be used to show that this
invariant distribution is unique.

######  Theorem 2.

[10] If a Markov process is Harris irreducible and has an invariant
probability $\mu^{*}$, then $\mu^{*}$ is the unique invariant probability and,
for any probability measure $\mu$ on S:

$\lim_{n\rightarrow\infty}\sup_{A\in\mathcal{B}(S)}\left|\frac{1}{n}\sum_{m-1}^{n}P^{m}\mu(A)-\mu^{*}(A)\right|=0$

Thus, if a Markov process is Harris irreducible, we know that any measure
converges to a unique invariant measure, $\mu^{*}$=, under the operator, $P$
“in the sense of Césaro”, i.e. the average of the sequence of measures
conerges to $\mu^{*}$. To get asymptotic stability of the invariant meausure
without this average, we need one more condition: aperiodicity.

######  Definition 6.

A Markov process is said to be strongly aperiodic if there is a set
$A_{0}\subset S$ with $\nu(A_{0})=1$ and an $\alpha>0$ such that for all $x\in
A_{0}$, $A$ a measurable subset of $S$, one has

$p(x,A)\geq\alpha\nu(A)$

######  Theorem 3 ([8], Theorem C9.2).

If a Markov process is Harris irreducible and strongly aperiodic with
invariant measure $\mu^{*}$ then, for any $\mu\in\mathcal{B}(S)$,

$\lim_{n\rightarrow\infty}\delta(P^{n}\mu,\mu^{*})=0$

where $\delta$ is the total variation distance

$\delta(P,Q)=\sup_{B\in\mathcal{B}}|P(B)-Q(B)|.$

As an example of the ideas presented in this section, let us consider a
version of the stochastic logistic map as examined in [11]. In this paper,
Bhattacharya and Rao consider the stochastic logistic map

$X_{n+1}=\lambda_{n+1}X_{n}(1-X_{n})$

where $\lambda$ is distributed according to a Bernoulli distribution, and
takes two possible values, $\alpha$ and $\beta$ (where $\alpha<\beta$). Thus,
the process switches randomly between two quadratic maps. Bhattacharya and Rao
demonstrate that if $0\leq\alpha<\beta\leq 1$, then $\delta_{0}$, the Dirac
measure supported on 0, is the unique invariant probability on [0,1]. If
$1<\alpha<\beta\leq 2$ then there is a unique invariant probability on (0,1),
and this probability is nonatomic. If $2<\alpha<\beta<1+\sqrt{5}$ and
$\alpha\in[8/\beta(4-\beta),\beta)$, then there exists a unique invariant
probability on (0,1) which is nonatomic and has its support contained in
$[1/2,\lambda/4]$.

##  3 Proof of the asymptotic stability of the invariant measure in the
chaotic regime.

In section 2, we outlined results that give conditions for the existence of an
invariant distribution of the stochastic logistic map, and for the stability
of this invariant distribution. In this section we use these results to prove
that, when parameters are supported in the chaotic regime, the random logistic
map has a unique asymptotically stable invariant measure.

Recall that the logistic family: $f_{\lambda}(x)=\lambda x(1-x)$ undergoes a
period-doubling cascade as $\lambda$ increases from 1 to approximately
3.56995. For higher values of $\lambda$, the family exhibits stable periodic
orbits of odd period, following the Sarkowski sequence, culminating with a
stable period-3 orbit when $\lambda=1+\sqrt{8}$. When $\lambda$ is greater
than this value the maps exhibit chaotic behavior on a Cantor set of parameter
values. In this section, we are interested in randomly perturbing the values
of $\lambda$ in this chaotic regime: $\lambda\in(1+\sqrt{8},4)$. The logistic
family is well-studied in the literature, and more details on its dynamics can
be found in, for example, [12]

In [13], Athreya and Dai consider the stochastic logistic map and show that if
$E(\log\lambda_{1})>0,E(\log(4-\lambda_{1})<\infty$, and if the process is
Harris irreducible, then the probability distribution converges in the sense
of Césaro to a unique probability distribution on (0,1). They also demonstrate
that the Markov process $\\{X_{n}\\}$ is Feller continuous, which guarantees
the existence of the invariant distribution (but not its stability) by Theorem
1. We shall consider the case when $\lambda\sim\mathcal{U}(3.87,4).$

######  Remark 1.

Let $\lambda_{i}\sim\mathcal{U}(3.87,4)$. Then we claim that
$E(\log\lambda_{1})>0$ and $E(\log(4-\lambda_{1}))<\infty$. Notice that if
$\lambda_{1}\sim\mathcal{U}(3.87,4)$, then $E(\lambda_{1})=3.935>0$.
Similarly,

$\displaystyle E(\log(4-\lambda))$ $\displaystyle=$
$\displaystyle\int_{3.87}^{4}\frac{1}{0.13}\lambda\log(4-\lambda_{1})d\lambda$
$\displaystyle\approx$ $\displaystyle-12,$

the last step computed numerically via Mathematica. Thus
$E(\log(4-\lambda_{1}))$ is finite.

This remark implies that if the process $\\{X_{n}\\}$ is Harris irreducible,
then there exists a unique probability distribution on $(0,1)$ to which other
probability distributions converge to in the sense of Cesaro.
The following results prove the Harris irreducibility and aperiodicity of the
Markov process, which then gives the asymptotic stability of the invariant
measure. Lemmas 1 and 2 build toward Theorem 4, which is a local
irreducibility result. Next, we present Theorem 5, which gives the Harris
irreducibility of the Markov process, and finally Theorem 6, which shows that
the process is strongly aperiodic and, hence, has a unique asymptotically
stable invariant measure.

######  Lemma 1.

Let $\lambda_{i}\sim\mathcal{U}(3.87,4)$ for all $i\geq 1$, $\alpha$ a fixed
element in $(3.87,4)$, and $x^{*}$ such that $x^{*}=\alpha x^{*}(1-x^{*})$.
(Equivalently, $x^{*}=1-\frac{1}{\alpha}$.) Then there exists $\eta>0$ such
that for all $x\in I=(x^{*}-\eta,x^{*}+\eta)$ and $A\subset I$ with $m(A)>0$,
then

$p(x,A)=\hbox{Prob}\\{X_{1}\in A\,|\,X_{0}=x\\}>0.$

That is, there is a neighborhood $I$ around the fixed point for which every
$x\in I$ has a positive probability of visiting each subset of $I$ of positive
measure.

###### Proof.

Fix $\alpha\in(3.87,4)$, and let $x^{*}\in(0,1)$ be such that $\alpha
x^{*}(1-x^{*})=x^{*}$. Notice that for any $\eta>0$, $x\in S$ and $A\subset
I\equiv(x^{*}-\eta,x^{*}+\eta)$,

$\sup\left\\{\frac{A}{x(1-x)}\right\\}\leq\frac{x^{*}+\eta}{\inf\\{x(1-x)|x\in
I\\}}\equiv a(\eta)$

and

$\inf\left\\{\frac{A}{x(1-x)}\right\\}>\frac{x^{*}+\eta}{\sup\\{x(1-x)|x\in
I\\}}\equiv b(\eta)$

Since $h(x)=x(1-x)$ is continuous at $x^{*}$ and $\alpha h(x^{*})=x^{*}$, it
follows that there exists $\eta>0$ such that $a(\eta)<4$ and $b(\eta)>3.87$.
Thus, $\frac{A}{x(1-x)}\subset(3.87,4)$, and since $A$ has positive measure,
so does $\frac{A}{x(1-x)}$, and therefore $P(x,A)>0$.

∎

######  Remark 2.

Note that while we take the distribution of $\lambda$ on $(3.87,4)$ to be
uniform, the same result will hold if we replace the uniform distribution with
a distribution that has non-vanishing absolutely continuous part.

######  Lemma 2.

Let $I$ be defined as in Lemma 1. For each $x\in(0,1)$, there exists
$n\in\mathbb{N}$, and a finite sequence
$\\{\alpha_{1},\alpha_{2},\ldots,\alpha_{n}\\}\subset(3.87,4)$ such that
$Y_{n}\in I$, $Y_{0}=x$, and $Y_{j+1}=f(Y_{j},\alpha_{j+1})$ for
$j=0,1,\ldots,n-1$. (That is, all $x\in(0,1)$ can get into $I$ within a finite
number of steps.)

###### Proof.

Fix $\alpha\in(3.87,4)$; recall then that the logistic map with parameter
$\alpha$ is _topologically transitive_ ; that is, there exists an orbit
$\\{y_{1},y_{2},y_{3},\ldots\\}$ that is dense in (0,1). Let
$0<\varepsilon<\min\\{\alpha-3.87,4-\alpha\\}$. We claim that, given
$x\in(0,1)$ there exists $\alpha^{*}\in(3.87,4)$ and $N_{1}\in\mathbb{N}$ such
that $f(\alpha^{*},x)=y_{N_{1}}$. In order for this to happen, we need

$\alpha^{*}=\frac{\alpha y_{N_{1}-1}(1-y_{N_{1}-1})}{x(1-x)}$

By the denseness of the orbit $\\{y_{i}\\}_{i=1}^{\infty}$, for all $\delta>0$
there is some $N_{\delta}\in\mathbb{N}$ such that $|x-y_{N_{\delta}}|<\delta.$
We can guarantee that $\alpha^{*}\in(3.87,4)$ by making an appropriate choice
for $\delta$, and letting $N_{1}=N_{\delta}$. Notice that

$\displaystyle|\alpha^{*}-\alpha|$ $\displaystyle=$
$\displaystyle\left|\frac{\alpha
y_{N_{1}-1}(1-y_{N_{1}-1})}{x(1-x)}-\alpha\right|$ $\displaystyle=$
$\displaystyle\alpha\left|\frac{y_{N_{1}-1}(1-y_{N_{1}-1})}{x(1-x)}-1\right|$
$\displaystyle<$
$\displaystyle\alpha\left|\frac{(x+\delta)(1-(x-\delta))}{x(1-x)}-1\right|$
$\displaystyle=$
$\displaystyle\alpha\left|\frac{x(1-x)}{x(1-x)}+\delta\left(\frac{1+\delta}{x(1-x)}\right)-1\right|$
$\displaystyle=$
$\displaystyle\alpha\delta\left|\frac{1+\delta}{x(1-x)}\right|$

which can thus be made arbitrarily small (and in particular, less than
$\varepsilon$) by taking $\delta$ sufficiently small.
Because $I$ has positive measure, the dense orbit $\\{y_{n}\\}_{n=1}^{\infty}$
must visit $I$ infinitely often. Therefore, there exists $N_{2}\geq N_{1}$
such that $y_{N_{2}}\in I$. Thus, the sequence
$\\{\alpha^{*},\alpha,\alpha,\ldots,\alpha\\}$ will map $x\in(0,1)$ into $I$
within a finite ($N_{2}-N_{1}+1$) number of iterations. ∎

######  Theorem 4.

$\\{X_{n}\\}$ is $m$-irreducible on $I$ where $m$ is Lebesgue measure.

###### Proof.

Fix $x\in(0,1)$ and let $A\subset I$ with $m(A)>0$. By continuity of $f$ and
by Lemma 2, for all $x\in(0,1)$, there exists $\varepsilon>0$ such that for
$\\{\alpha_{i}\\}_{i=1}^{n}$ as given in Lemma 2, if
$|\lambda_{i}-\alpha_{i}|<\varepsilon$ for $i=1,2,\ldots,n$, where $X_{0}=x$,
$X_{j+1}=f(X_{j},\lambda_{j+1})$ for $j=0,1,\ldots,n-1$, then $X_{n}\in I$.
Notice then that

$p(X_{n}\in
I|X_{0}=x)\geq\prod_{i=1}^{n}p(|\lambda_{i}-\alpha_{i}|<\varepsilon)>0.$

By Lemma 1, $p(X_{n+1}\in A|X_{n}\in I)>0.$ Thus,

$p(X_{n+1}\in A|X_{0}=x)=p(X_{n+1}\in A|X_{n}\in I)\cdot p(X_{n}\in
I|X_{0}=0)>0.$

Therefore $\\{X_{n}\\}$ is $m$-irreducible on $I$. ∎

######  Remark 3.

We can equivalently think of $\\{X_{n}\\}$ as Harris irreducible on $S$ with
reference measure $\phi(\cdot)=m(\cdot\cap I)$, where $m$ is Lebesgue measure.

######  Theorem 5.

As above, let $S=[0,1]$.

$X_{n+1}=\lambda_{n+1}X_{n}(1-X_{n}),\,\,\,n\geq 0$

with $\\{\lambda_{n}\\}_{n\geq 1}$ being i.i.d.,
$\lambda_{i}\sim\mathcal{U}(3.87,4)$, and $X_{0}$ an independent $S$-valued
r.v. Then $\\{X_{n}\\}$ is irreducible with respect to the Lebesgue measure
$m$ on $S$.

###### Proof.

By a result of Graczyk and Swiatek [14], the set of all $\lambda$ in $(0,4)$
such that $f(\cdot,\lambda)$ has a stable periodic orbit is dense. Thus, there
is a $\alpha\in(3.87,4)$ and $0<p<1$ such that $x^{*}$ is a stable periodic
point for $f(x,\alpha)=\alpha x(1-x)$. Let $\delta$ be such that
$(\alpha-\delta,\alpha+\delta)\subset(3.87,4)$.

Let $K=\\{x\in(0,1)\,\,|\,\,\omega_{\alpha}(x)=\gamma(x^{*})\\}$, where
$\omega_{\alpha}$ is the limit set of the orbit of $x$ under $f(\cdot,\alpha)$
and $\gamma(x^{*})$ is the orbit of the periodic point $x^{*}$. In [15]
(Chapter V, Theorem 1.3) de Melo and van Strien have shown that $m(K)=1$ for
unimodal maps, where $m$ is Lebesgue measure. Thus let
$I=(x^{*}-\eta,x^{*}+\eta)$ as given by Theorem 1. Then, for all $x\in K$,
$f^{(n)}(x,\alpha)\in I$ for some large $n$. Thus, by continuity of $f$,

$p_{x}(X_{n}\in I\text{ for some }n\geq 1)>0$ (3)

for all $x\in K$.

We will next show that for any $x\in(0,1)$, $P_{x}(X_{1}\in K)>0$. $x\in(0,1)$
implies $\alpha x(1-x)\in(0,1)$, and since $m(K)=1$,

$\displaystyle m\left(\frac{K}{h(x)}\cap(3.87,4)\right)$ $\displaystyle\geq
m([0,\alpha]\cap(3.87,4))$ $\displaystyle\geq\delta>0.$

Thus,

$p\left(\lambda_{1}\in\frac{K}{h(x)}\cap(3.87,4)\right)\geq\frac{\delta}{0.13}>0$

and

$p_{x}(X_{1}\in K)=p\left(\lambda_{1}\in\frac{K}{h(x)}\right)\geq
p\left(\lambda_{1}\in\frac{K}{h(x)}\cap(3.87,4)\right)>0,$

establishing irreducibility. ∎

We now have that, since $E(\lambda)>0$ and $E(\log(4-\lambda))<\infty$, and
the process is Harris irreducible, by [13], the process $\\{X_{i}\\}$ has a
unique invariant measure that other measures converge to in the sense of
Cesaro. We conclude this section by demonstrating that the process is strongly
aperiodic, and thus the invariant distribution is asymptotically stable with
respect to the total variation norm.

######  Theorem 6.

Let $S=[0,1]$.

$X_{n+1}=\lambda_{n+1}X_{n}(1-X_{n}),\,\,\,n\geq 0$

with $\\{\lambda_{n}\\}_{n\geq 1}$ being i.i.d.,
$\lambda_{i}\sim\mathcal{U}(3.87,4)$, and $X_{0}$ an independent $S$-valued
r.v. Then $\\{X_{n}\\}$ is strongly aperiodic, and thus has an asymptotically
stable invariant distribution.

###### Proof.

Let $\lambda_{min}=3.87$, and let
$A_{0}=(1-\frac{1}{\lambda_{min}},1-\frac{1}{4})$ (that is, $A_{0}$ is the set
of fixed points for the deterministic systems corresponding to the possible
parameter $\lambda$ values of the logistic map). Let $L$ be the length of
$A_{0}$ (that is, $L_{1}=(1-\frac{1}{4})-(1-(\frac{1}{\lambda}_{min})$ and let
$L_{2}$ be the length of interval $(\lambda_{min},4).$ Let $\nu$ be normalized
Lebesgue measure on $A_{0}$. Then $A_{0}$, by definition, has measure 1, and
the only sets of positive measure are those that overlap $A_{0}$ in a set of
positive Lebesgue measure. Thus, notice that

$\nu(A)=\frac{1}{L_{1}}\int_{y\in A}dy$

Fix $x\in A_{0}.$ Then

$\displaystyle p(x,A)$ $\displaystyle=$ $\displaystyle
p\left(\lambda\in\frac{A}{x(1-x)}\right)$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{L_{2}}\int_{\frac{A}{x(1-x)}}d\lambda$
$\displaystyle=$ $\displaystyle\frac{1}{L_{2}}\int_{A}\frac{1}{x(1-x})dy$
$\displaystyle=$ $\displaystyle\frac{1}{L_{2}}\frac{1}{x(1-x)}\int_{A}dy$
$\displaystyle=$ $\displaystyle\frac{1}{L_{2}}\cdot\frac{1}{x(1-x)}\cdot
L_{1}\nu(A)$ $\displaystyle\geq$ $\displaystyle
4\cdot\frac{L_{1}}{L_{2}}\nu(A),$

since $x(1-x)\leq\frac{1}{4}$ for all $x\in[0,1]$. Thus, letting
$\alpha=4\cdot\frac{L_{1}}{L_{2}}$ gives us the desired result. Therefore, the
process $\\{X_{n}\\}$ is strongly aperiodic, and has an asymptotically stable
invariant distribution. ∎

##  4 Simulations and Discussion

Suppose $X$ is a random variable taking values in $[0,1]$. If $Pr\\{X\leq
x\\}=Pr\\{f(X)\leq x\\}$ then we say that the distribution of $X$ is invariant
under $f$. In this way, we can discuss invariant distributions for the
deterministic logistic map, by looking at distribtuions invariant under
$f_{\lambda}(x)=\lambda x(1-x)$. In [16], Ulam and von Neumann proved that the
continuous invariant distribution for the deterministic logistic map with
$\lambda=4$ is a beta distribution with $\alpha=\beta=0.5$. This function is
symmetric on the interval (0,1), and goes to infinity at the endpoints of its
support. Thus, the invariant distribution when $\lambda=4$ is given by the
function

$f^{*}(x)=\begin{cases}\frac{1}{\pi\sqrt{x(1-x)}}&0<x<1\\\ 0&\hbox{
else}\end{cases}$ (4)

See Figure 1 for a comparison of this beta density and an empirical
distribution.

Figure 1: Histogram of the beta distribution with $\alpha=\beta=0.5$. The
solid line shows the theoretical distribution given in Equation (4); the bars
are a histogram generated by sampling 100,000 initial values uniformly on
$(0,1)$

However, invariant densities for other values of $\lambda$ are not, in
general, symmetric. In [17], Hall and Wolff provided numerical estimations for
the invariant distributions for the deterministic logistic map for various
values of $\lambda$. They demonstrated that any continuous invariant density
will have poles (singularities) at the boundary of its support, as well as
additional poles at the interior of the support in the case $\lambda<4$. They
also show that these invariant distributions were not symmetric. Figure 2
shows numerical approximations to theses deterministic invariant densities for
$\lambda=3.87$ and $\lambda=3.95$ These numerical approximations are
histograms of 100,000 sampled initial values iterated 2,000 times. They
indicate the singularities in the theoretical densities described in Remark 5
of [17]. These singularities, or poles, correspond to the orbit of the
critical point, $f_{\lambda}(1/2)$. In addition, Hall and Wolff observed that
the invariant density will be supported on the interval
$(\frac{1}{4}\lambda^{2}(1-\frac{\lambda}{4}),\frac{1}{4}\lambda)$. This is
due to the fact that any point in the interval $(\frac{1}{4}\lambda,1)$ will
be mapped into the interval $(0,\frac{1}{4}\lambda^{2}(1-\frac{\lambda}{4}))$,
and any point in $(0,\frac{1}{4}\lambda^{2}(1-\frac{\lambda}{4})$ will be
mapped to a point in
$[\frac{1}{4}\lambda^{2}(1-\frac{\lambda}{4}),\frac{1}{4}\lambda]$. Therefore,
after finitely many iterations, any point $x\in(0,1)$ must be mapped into
$[\frac{1}{4}\lambda^{2}(1-\frac{\lambda}{4}),\frac{1}{4}\lambda]$, and so the
invariant density must be supported on this subinterval. Note that as
$\lambda\rightarrow 4$,
$\frac{1}{4}\lambda^{2}(1-\frac{\lambda}{4})\rightarrow 0$ and
$\frac{1}{4}\lambda\rightarrow 1$, so the interval of support widens as
$\lambda$ approaches 4.

Figure 2: Histogram of the simulated invariant distribution for the
deterministic logistic map with $\lambda=3.87$ (left) and $\lambda=3.935$
(right) for 100,000 sampled values.

To numerically explore convergence to the invariant distribution in the
stochastic case, we ran several simulations with 100,000 initial values
sampled randomly according to a variety of distributions:

  1. 1.

A discrete distribution supported on the set $\\{0.11,0.33,0.55,0.6,0.78\\}$
with probabilities $\\{0.196,0.140,0.233,0.322,0.107\\}$ respectively (Figure
3, left).

  2. 2.

A discrete distribution supported on $\\{0.25,0.5,0.75\\}$ each with
probability 1/3.

  3. 3.

A truncated exponential distribution on (0,1) with parameter $r=1.2$ (Figure
3, right).

  4. 4.

A truncated gamma distribution on (0,1) with shape parameter 3 and scale
parameter 1.

  5. 5.

A truncated normal distribution on (0,1) with mean 0.5 and standard deviation
0.3.

  6. 6.

A truncated $t$-distribution on (0,1) with 1 degree of freedom.

  7. 7.

A uniform distribution on (0,1) (Figure 3, center).

In each case, after 2,000 iterations, we obtained a histogram very similar to
that presented in Figure 4: a bi-modal, asymmetric distribution that appears
to be supported on the entirety of the unit interval, without singularities
(poles) in the interior. This is noteworthy, particularly because we simulated
using both continuous and discrete starting distributions. We examined
distributions that were left-skewed, right-skewed, and symmetric. Each
simulated initial distribution appears to converge to the same distribution
after 2,000 iterations, aligning with the global asymptotic stability of the
invariant distribution established theoretically earlier in this paper. We
also note that, while the distribution demonstrates peaks toward the sides of
the interval, it does not exhibit the poles at the boundary of the interval in
the same manner that the invariant distributions for the deterministic maps
do. However, this seems reasonable, as the interval of support of the
invariant distribution is dependent on $\lambda$, as discussed above, and
relatively fewer values of $\lambda$ will map points into ranges closer to 0
and 1, so it makes sense that the distribution will decrease near the boundary
of $[0,1]$

Figure 3: Starting distributions of 100,000 starting values sampled on the set
$\\{0.11,0.33,0.55,0.6,0.78\\}$ with probabilities
$\\{0.196,0.140,0.233,0.322,0.107\\}$ respectively (left), the uniform
distribution on (0,1) (middle), and the truncated exponential distribution
with parameter $r=1.25$ (right). Figure 4: Empirical distribution of
$\\{X_{n}\\}$ after $n=2,000$ iterations, where 100,000 initial values were
sampled from an exponential distribution with parameter $r=1.25$. Similar
results were obtained using other initial distributions (results not shown).

##  5 Conclusion

In conclusion, we have illustrated that the stochastic logistic map,
characterized by parameter values uniformly distributed in the chaotic regime,
exhibits Harris irreducibility. Consequently, the previously established
invariant distribution stands as the unique invariant distribution. This
result also implies that all initial distributions converge to the invariant
distribution in the Cesàro sense. Moreover, we have demonstrated the strong
aperiodicity of the stochastic logistic map, indicating that initial
distributions converge to the invariant distribution with respect to the total
variation distance—a more robust and rigorous result. We conducted simulations
of the stochastic logistic map to visualize the invariant distribution.
Remarkably, the empirically derived invariant distribution exhibits
qualitative distinctions when compared to the invariant distributions of the
deterministic logistic map for values of $\lambda$ in the chaotic regime. In
particular, the simulated distribution does not appear to have a pole at
$x=0$, which contrasts with the work in [17], which demonstrates that the
deterministic invariant density will have a pole at the boundary of its
support.
It should be noted that this work does not seem to depend on the $\lambda$
values being uniformly distributed on (3.87,4) and indeed, it appears that the
results will hold for any absolutely continuous distribution on (3.87,4).

In future work, we plan to provide a theoretical basis for these empirical
observations. In particular, we hope to better describe the asymptotically
stable invariant distribution in the case that the $\lambda$-values are
uniformly distributed (the case illustrated in Figure 4). We also plan to
provide similar theoretical results when the parameter has other
distributions, such as an exponential or discrete distribution.

## References

  *   * May [1976] May, R.M.: Simple mathematical models with very complicated dynamics. Nature 261, 459–467 (1976)
  * Haskell, C. and Sacker [2005] Haskell, C., Sacker, R.J.: The stochastic beverton-holt equation and the m. neubert conjecture. Journal of Dynamics and Differential Equations 17(4), 825–844 (2005)
  * Descheemaeker and de Buyl [2020] Descheemaeker, L., Buyl, S.: Stochastic logistic models reproduce experimental time series of microbial communities. eLife 9, 55650 (2020) https://doi.org/10.7554/eLife.55650
  * Satoh and Yamada [2001] Satoh, D., Yamada, S.: Discrete equations and software reliability growth models. In: Proceedings 12th International Symposium on Software Reliability Engineering, pp. 176–184 (2001). https://doi.org/10.1109/ISSRE.2001.989471
  * Nishimura et al. [2004] Nishimura, K., Stachurski, J., et al.: Discrete time models in economic theory (2004)
  * Bhattacharya and Majumdar [2004] Bhattacharya, R., Majumdar, M.: Stability in distribution of randomly perturbed quadratic maps as Markov processes. The Annals of Applied Probability 14(4), 1802–1809 (2004) https://doi.org/10.1214/105051604000000918
  * Kifer [1986] Kifer, Y.: Ergodic Theory of Random Transformations. Birkhauser, Boston, Basel, Stuttgart (1986)
  * Bhattacharya and Majumdar [2007] Bhattacharya, R., Majumdar, M.: Random Dynamical Systems. Cambridge University Press, Cambridge (2007)
  * Lasota, A. and Mackey [1994] Lasota, A., Mackey, M.C.: Chaos, Fractals, and Noise, 2nd edn. Springer, Berlin, Heidelberg, New York (1994)
  * Nummelin [1984] Nummelin, E.: Miscellaneous limit theorems for Harris recurrent Markov chains. Cambridge Tracts in Mathematics, pp. 126–140. Cambridge University Press, ??? (1984). https://doi.org/10.1017/CBO9780511526237.008
  * Bhattacharya and Rao [1993] Bhattacharya, R.N., Rao, B.V.: Random Iterations of Two Quadratic Maps, pp. 13–22. Springer, New York, NY (1993). https://doi.org/10.1007/978-1-4615-7909-0_3 . https://doi.org/10.1007/978-1-4615-7909-0_3
  * Holmgren [1996] Holmgren, R.: A First Course in Discrete Dynamical Systems, 2nd edn. Springer, Berlin, Heidelberg, New York (1996)
  * Athreya and Dai [2000] Athreya, K.B., Dai, J.: Random logistic maps. i. Journal of Theoretical Probability 13, 595–608 (2000)
  * Graczyk and Swiatek [1997] Graczyk, J., Swiatek, G.: Generic hyperbolicity in the logistic family. Annals of Mathematics 146(1), 1–52 (1997). Accessed 2023-10-07
  * de Melo and van Strien [1993] Melo, W., Strien, S.: One-dimensional Dynamics. Springer, ??? (1993)
  * Ulam and von Neumann [1947] Ulam, S.M., Neumann, J.: On combination of stochastic and deterministic processes. Bulletin of the American Mathematical Society 53, 1120 (1947)
  * Hall and Wolff [1995] Hall, P., Wolff, R.C.L.: Properties of invariant distributions and lyapunov exponents for chaotic logistic maps. Journal of the Royal Statistical Society. Series B (Methodological) 57(2), 439–452 (1995). Accessed 2023-11-11


# Multi-Genre Music Transformer - Composing Full Length Musical Piece

Abhinav Kaushal Keshari (Purdue University)

###### Abstract

In the task of generating music, the art factor plays a big role and is a
great challenge for AI. Previous work involving adversarial training (Dong et
al., 2018) to produce new music pieces and modeling the compatibility (Huang
et al., 2021) of variety in music (beats, tempo, musical stems) demonstrated
great examples of learning this task. Though this was limited to generating
mashups or learning features from tempo and key distributions to produce
similar patterns. Compound Word Transformer (Hsiao et al., 2021) was able to
represent music generation task as a sequence generation challenge involving
musical events defined by compound words. These musical events give a more
accurate description of notes progression, chord change, harmony and the art
factor. The objective of the project is to implement a Multi-Genre Transformer
which learns to produce music pieces through more adaptive learning process
involving more challenging task where genres or form of the composition is
also considered. We built a multi-genre compound word dataset, implemented a
linear transformer (Katharopoulos et al., 2020) which was trained on this
dataset. We call this Multi-Genre Transformer, which was able to generate full
length new musical pieces which is diverse and comparable to original tracks.
The model trains 2-5 times faster than other models discussed.

Artificial Intelligence

##  1 Related Work

Despite achieving great success in generation challenges using Artificial
Intelligence in Natural Language Generation (NLG) there is a factor of art
that still makes them different from human like performance. In terms of NLG
we can relate it to something like the difference between computer generated
article and a piece of art like novels, biography, etc. For music art factor
always come into account and despite able to produce musical compositions
through Adversarial networks or mixing stems using supervised learning the
solution still is very different from an original piece of music which we
discuss below.

###  1.1 Music Generation using GANs

Generative adversarial networks (GANs) have provided significant progress in
producing text, videos and images. Similar efforts have been made to bring
neural networks to artistic domain of music. MuseGAN(Dong et al., 2018)
brought a novel model for generating multi-track music.

Until 2018, the progress in using AI to compose music had been able to produce

  * •

Single-track (monophonic) music

  * •

Multi-track (polyphonic) music by combining several monophonic melodies in
chronological order

Music usually being an art involving multiple instruments played together
requires music to be multi-track and because music notes are made up of
chords, arpeggios or melodies the idea of using a chronological order setting
prevents it from being generalized.

The paper(Dong et al., 2018) address this challenge in generalising real music
by discussing current technical lacks in neural network models and how it
relates to the real world music.

  1. 1.

Music is an art of time and has characteristics of coherence, rhythm, tension
and emotion flow. This requires it to have a Temporal Model.

  2. 2.

Music compositions usually involves different instruments interacting with one
another making the compositions to be harmonic. To solve this issue a Composer
Model is required.

  3. 3.

Musical notes are built of chords, arpeggios or melodies and how they unfold
over time; thus introducing chronological generation of notes is not suitable.
To address this the paper introduces using bars (segment of time) instead of
notes as the basic unit for composition. And then generate music bar by bar
using transposed convolutional neural networks to learn translation-invariant
patterns.

The paper(Dong et al., 2018) makes contributions in terms of both ability to
artificially compose realistic music and use of generative adversarial
framework with temporal and composition models. In short the contributions
are:

  * •

First GAN based model for generating multi-track sequence.

  * •

First model which can generate multi-track polyphonic music.

  * •

Same model can be used as a music accompaniment.

  * •

Creates a new Lakh Pianoroll Dataset (LPD) for multitrack piano-rolls

  * •

For future work metrics in the domain of artificial music a new set of
objective metrics are proposed.

MuseGAN model proposed considers two sub-network generator Gtemp (temporal
structure generator) and Gbar (bar generator) making the overall generator:

$G(z)=\left\\{G_{bar}(G_{temp}(z)^{(t)})\right\\}^{T}_{t=1}$

where z is the input noise vector. The strength of the model is the ability to
generate samples having chord like intervals (learning features from temporal
model) and melodies involving pitch overlap among guitar, piano and strings
(learning features from composer model).

The model introduces multi-track by modeling interdependency of tracks by
proposing 3 different generator model (Jamming, Composer and Hybrid), but the
author brings up these based on the understanding of pop music composition.
This possibly restricts the generator to explore on a broad spectrum of music
and prevents it from being generalised. Also worth mentioning is that the work
relies on multi-track interdependency, but misses to study about the
compatibility of these tracks which can significantly increase the quality of
music being generated. We will see this issue being addressed in the next
paper.

###  1.2 Modeling the Compatibility of Stem Tracks to Generate Music
Mashups(Huang et al., 2021)

Source separation(Jansson et al., 2017; Défossez et al., 2019) makes it
possible to generate a music mashup with isolated stems like vocals, drums,
piano, etc. The challenge lies in producing music which has compatibility
between these stems. This paper creates a mashup generation pipeline and
trains a model to predict the compatibility by automatically learning to
adjust key and tempo (characteristics of quality mashups in real world).

General models trained for harmonic compatibility (Bernardes et al., 2017;
Macas et al., 2018) fails to consider subtle features or surprise mixes of
disparate samples which is quite common in this art domain. Other issue that
arises is audio compatibility models like Neural Loop Combiner (Chen et al.,
2020) having lack of vocal source and variety of genres.

The authors designed a self supervised learning model by recombining the
original combination of stems before source separation to serve as examples of
ground truth. To avoid highly polarized model, semi-supervised learning was
introduced which included producing several random mashups by mixing different
stems and treated them as unlabeled instances. Label smoothing regularization
for outliers (Zheng et al., 2017) was used to assign uniform distribution to
the unlabeled data for loss computation. This helps in regularization.

The final architecture consists of 3 modules:

  1. 1.

Music Source Separation: Uses MSS algorithm (Jansson et al., 2017) to get
different stems vocals, drums, bass and other.

  2. 2.

Mashup Database (MashupDB): Using Madmom (Böck et al., 2016) different
features from the music clips are extracted like key, tempo and downbeat
information. Using these features and separate stem combinations a mashup
database is created which will act as either harmonic or percussion stem
candidates for mashup generation process.

  3. 3.

Mashup Generation: It uses candidate stems from MashupDB and adjusts key and
tempo to produce mashups within 3 conditions - original, matched and
unmatched.

The model (Huang et al., 2021) is defined by $p(y|V,H,P)$ where $V$, $H$, and
$P$ are input signals for respective stems vocal, harmonic, and percussion.
The output probability p is used as the mashup compatibility and
$y\in\left\\{0,1\right\\}$ stating good or bad.

The model (Huang et al., 2021) implementation tries to mimic learning
compatibility for producing new mashups and provides objective and subjective
evaluation by cross validation among multiple different datasets. This
technique becomes easier because of the ability of the model to extract
different stems and features and build its own mashup candidates. This also
makes the model training process not dependent on human labeled data. The
model is also robust as negative data is added along with positive data for
supervised learning. The range of music coverage is also extensive and the
source separation step makes it easier for the model to be extended to
different genres for training.

But the current model design lacks the effective embedding of different stems
while producing a mashup and makes it dependent on tuning of key and tempo.
Currently the implementation comes up with fixed range of key and tempo
difference for compatibility and does not explain in detail how they came up
with these numbers. Although defining a range prevents large pitch shifting
and time stretching. Additionally the results of the model ranks positive
labeled data (original) over unlabeled data which might lead to concerns of
flexibility. Another major challenge of the model is the large training time
which is around 3 days using an NVIDIA Tesla-V100 GPU whereas using
transformer model significantly reduces the training time.

###  1.3 Music Transformers

With state-of-the art neural network we managed to learn features in music by
defining certain rules on matching tempo, beats or compatibility. In the
previous paper we also tried to learn compatibility with the help of
supervised learning. The model though suffered with bias as compatibility was
favoured for matched key or tempo and also lacks generalization. Compound Word
Transformer (Hsiao et al., 2021) considers music as sequence of events and
uses a Transformer (neural sequence model) (Vaswani et al., 2017) to generate
a new musical sequence.

A musical note can be described by note’s pitch, chord, bar, duration,
velocity (dynamics), placement (onset time). If we consider these as tokens we
can then define music as sequence of tokens and these tokens are a part of
pre-defined vocabulary. As music is multi-faceted a particular type of token
can capture only a certain feature like melody, rhythm, harmony. All the
neural networks until now treated these tokens as equal and thus lacked
heterogeneity.

Compound Word Transformer (Hsiao et al., 2021) generates music in a
conceptually different way as it allows tokens to be of specific types and let
them have their own properties. Tokens can be of note type (pitch, duration)
or metric type (beginning of new beat, bar). We then defines a musical event
by combination of such tokens which allows to capture co-occurrence
relationship among the tokens. This combination of tokens are termed as
compound words. So, now we can represent a music piece $(X)$ as a sequence
$(S)$ of compound words $(cp)$ or $S=g(X)=\left\\{cp_{t}\right\\}^{T}_{t=1}$
where $g(.)$ is the conversion function to convert music into time-ordered
sequence of musical events and $T$ is the length of the music sequence.

Figure 1: Graphical Representation of Music Space

Theoretically, the model learns over discrete-time dynamic directed
hypergraphs. Consider a graph $G=(V,E)$ (Figure 1) the vertices $(V)$ are
tokens and edges $(E)$ are sequence of token. Collection of vertices can be
defined as a compound word and hyperedge in this graph represents sequence of
compound words. In figure 1 $v1,v2,v5$ are the tokens and the edge $E1$
defines a sequence of tokens whereas $e1,e2$ defines a hyperedge (connecting
more than 2 nodes). And transitioning from one hyperedge to another defines
the sequence of composition words which we are trying to learn.

Using a transformer we are trying to learn the next musical event or compound
word (combination of tokens). The self attention part of the transformer
learns the dependency among the elements in musical sequence and different
feedforward head is used for tokens of different type. In short the
implementation groups tokens to form compound words and then perform sequence
modeling in this sequence of compound words, the major contributions are:

  * •

Compose pop-piano music of full song length.

  * •

Compound word sequencing with linear transformer providing state-of-the-art
results in terms of quality with 5-10x faster training and inference time.

  * •

Music defined as Dynamic Directed Hypergraph.

Generating a new musical event or a group of tokens to be combined as a
compound word at each time step is the backbone of this model, but it relies
on assuming that no two musical events can occur together. The new hyperedge
generated by the Transformer decoder marks other tokens as [ignore] once an
event of a particular token type is detected. Can this limit the music
generation task? Additionally the model is trained using only pop music which
limits the expressing power of the transformer.

##  2 Implementation

Compound Word Transformer (Hsiao et al., 2021) was able to represent music
generation task as a sequence generation challenge involving musical events
defined by compound words. Leveraging this representation we implement a
neural model which learns to produce music pieces through more adaptive
learning process involving more challenging task where genres or form of the
composition is also considered. This adds the richness of music art in the
learning process of attention driven sequential learning. We will call this
model Multi-Genre Music Transformer and following are the steps involved for
implementing this:

  * •

Building Dataset: This involves generating compound word dictionary for songs
of different genres.

  * •

Implementing Transformer Model: We implement our Transformer class, the
training steps and the generation logic for inference.

  * •

Adaptive Learning: We allow our tuned model to be adaptable by training on a
smaller and multi-genre dataset.

###  2.1 Building Dataset

To be able to provide a more generalised learning process for our transformer
it needs to be trained with a piano roll dataset involving musical pieces of
variety of genres/style. The dataset should be based on compound words (Hsiao
et al., 2021) to represent different musical tokens as a combined unit for
sequence modeling which is different from traditional musical dataset (MIDI,
REMI).

Figure 2: Dataset Building Pipeline

This required us to build a dataset by selecting music clippings and
converting them to piano roll using Onsets and Frames (Hawthorne et al.,
2017). Extracting downbeat and beat information from these songs using madmom,
a music signal processing library (Böck et al., 2016). Finally representing
these metadata into a compound word representation using the dataset
generation scripts provided in the compound word transformer
repository111https://github.com/YatingMusic/compound-word-
transformer/blob/main/dataset/Dataset.md. This also adds on to the AILabs.tw
Pop1K7 dataset (Hsiao et al., 2021) which currently only includes pop music.
Figure 2 demonstrates the pipeline for creating a new dataset.

Following the pipeline above we managed to create a Compound Word (Hsiao et
al., 2021) dataset which involved piano roll for 150 musical pieces from 3
different genres including Electronic Dance Music (EDM), Indie and Hip-Hop.

###  2.2 Implementing Transformer Model

We implement a linear transformer(Katharopoulos et al., 2020) to address long
sequence dependency which is a very relevant factor in music generation due to
the presence of a context or a rhythm in the entire musical piece. Having an
independent feed-forward head in the Transformer Decoder allows to improve the
loss of independent tokens. This allows the model to scale for additional
perspective (like genre, form or involving a particular chord progression) in
the music by adding an additional token type. We implement our transformer
model in a generic way which allows user to define its own token sampling
model, token embedding model and these can be scalable for any number of token
types. The loss observed at each feed-forward head is shown in Figure 6. This
shows adding a new token (for genre/style/form) for model to learn can be
simply achieved by adding an independent feed-forward head for the same.

####  2.2.1 Token Embedding

Figure 3: Demonstrates how each token undergoes independent embedding before
combining with Positional Encoding. Here $T_{1},T_{2}...T_{k}$ are $K$
different tokens for our Transformer each having its own embedding function
and dimension. We are assuming the Transformer supports $K$ type of tokens.

The input to a transformer requires positional encoding added to the embedding
vector of our input sequence elements. As each element in our sequence is a
compound word (Hsiao et al., 2021) which is combined of different tokens, we
embed each token separately (allowing to have adaptive size) and then
concatenate them. Having an adaptive token size allows to use smaller
embedding dimension for a token type with smaller vocabulary and when we
concatenate all of these we get an embedding dimension of 512 for our model.
Refer to Figure 3 for detailed steps of token embedding.

####  2.2.2 Token Sampling

For inference, sampling plays a crucial role to avoid degeneration and improve
diversity. To avoid degeneration we follow Nucleus Sampling (Holtzman et al.,
2019), which is a stochastic temperature controlled process. This method
samples from the smallest subset of tokens whose cumulative probability mass
exceeds a threshold. We also had each token to have a separate sampling policy
by defining different threshold $p$ and different temperature parameter $\tau$
(Ackley et al., 1985) for reshaping the probability before sampling. We reused
the inference implementation from Compound Word Transformer (Hsiao et al.,
2021) and tweaked $\tau$ to have higher values for chord to allow more diverse
chord progressions. Figure 4 shows the sampling process and individual feed-
forward layer for each token in the transformer.

Figure 4: Transformer with N self-attention layers and independent feed-
forward head for each token. We first predict the Token Type for the
particular time-step and then perform a nucleus sampling before predicting the
remaining tokens.

###  2.3 Adaptive Learning

After defining the model, the next important step is to implement the training
steps. To support scalable token definition in our generalised transformer we
make the training steps modular and general to variable number of token types.
This allows easy addition of a new token and independently monitor gradient
descent optimization for the respective loss.

We trained our model in parallel for 2 different conditions. The first set of
training was performed on the original AILabs.tw Pop1K7 dataset (Hsiao et al.,
2021). The second set of training took into consideration to provide multi-
genre learning environment for the transformer as it involved training on a
dictionary that was generated from 3 different genres (EDM, Indie, Hip-Hop).

##  3 Evaluation and Results

To train a multi-genre transformer the primary objective was to provide it
with a dataset which is richer in variety than the original pop only dataset.
With the help of dataset building pipeline we managed to create a token set
which has a higher variance allowing the model to have a broader expressive
power. Figure 5 shows the comparison of tokens between the 2 datasets used.

Figure 5: Left image shows token distributions for the songs in the generated
multi-genre dataset and the right image shows similar distribution for
AILabs.tw Pop1K7 dataset (Hsiao et al., 2021).

After training the model for both the datasets we also observe (refer to
Figure 6) the individual token loss and total average loss is similar and
indicates the model converging. Additionally, the gradient descent is more
gradual using the multi-genre dataset displaying a more settled progression.

We trained the model with 12 self-attentions layers, 8 feed-forward heads with
model dimension of 512 and batch size of 4 for 180 epochs which took around
17hrs. Then using the trained model we generated 20 new full length musical
pieces with an average inference time of 12.56sec/song which is faster than
the compound-word transformer though having slightly less number of average
tokens per song. Table 1 shows a more detailed comparison.

Figure 6: Loss vs Epoch for different token types. The last plot corresponds
to the average loss for all different token types. Table 1: Quantitative
evaluation results for Multi-Genre Transformer and Compound Word Transformer.
Results for Compound Word Transformer comes from the implementation in the
paper (Hsiao et al., 2021).
Model | Training Time | GPU | Inference Time (/song) | Avg Tokens (/song)
---|---|---|---|---
Multi-Genre Transformer | 17 hrs | 9.8GB | 12.56 sec | 9190
Compound Transformer | 1.3 days | 9.5GB | 19.8 sec | 9546

For a qualitative evaluation of the musical pieces that were produced we
compare (Figure 7) the piano rolls of these with the piano rolls of original
tracks that were used to train the model.

[3pt]Original Songs
[3pt]Generated Songs

Figure 7: Piano roll of original and generated songs. We can see a rich and
complete content for the generated songs similar to some original tracks.

##  4 Conclusion

In this project we produce music as a sequence of musical events produced by a
trained Transformer. We leverage the definition of Compound Word (Hsiao et
al., 2021) to define musical event by grouping multiple tokens. This grouping
greatly reduces the size of our sequence and boosts long-range learning. This
also reduces the training and inference time for our model remarkably. We also
exploit the feature of each token having its independent feed-forward head for
prediction to make the model scalable for new token types that can be
introduced in our dictionary. This allows to add any new token for this
transformer very easily which can be used for musical form, chord progression,
etc. Additionally, we created an entire new dataset consisting of multi-genre
compound word dictionary and trained our model with this to provide it a more
adaptive learning environment. The compositions that were generated were
highly rich in musical events and were of good quality.

## References

  * Ackley et al. (1985) Ackley, D. H., Hinton, G. E., and Sejnowski, T. J.  A learning algorithm for boltzmann machines.  _Cognitive science_ , 9(1):147–169, 1985.  URL https://www.sciencedirect.com/science/article/abs/pii/S0364021385800124.
  * Bernardes et al. (2017) Bernardes, G., Davies, M. E., and Guedes, C.  A hierarchical harmonic mixing method.  In _International Symposium on Computer Music Multidisciplinary Research_ , pp. 151–170. Springer, 2017.  URL https://link.springer.com/chapter/10.1007/978-3-030-01692-0_11.
  * Böck et al. (2016) Böck, S., Korzeniowski, F., Schlüter, J., Krebs, F., and Widmer, G.  Madmom: A new python audio and music signal processing library.  In _Proceedings of the 24th ACM international conference on Multimedia_ , pp. 1174–1178, 2016.  URL https://dl.acm.org/doi/abs/10.1145/2964284.2973795.
  * Chen et al. (2020) Chen, B.-Y., Smith, J. B., and Yang, Y.-H.  Neural loop combiner: Neural network models for assessing the compatibility of loops.  _arXiv preprint arXiv:2008.02011_ , 2020.  URL https://arxiv.org/abs/2008.02011.
  * Défossez et al. (2019) Défossez, A., Usunier, N., Bottou, L., and Bach, F.  Music source separation in the waveform domain.  _arXiv preprint arXiv:1911.13254_ , 2019.  URL https://arxiv.org/abs/1911.13254.
  * Dong et al. (2018) Dong, H.-W., Hsiao, W.-Y., Yang, L.-C., and Yang, Y.-H.  Musegan: Multi-track sequential generative adversarial networks for symbolic music generation and accompaniment.  In _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018\.  URL https://www.aaai.org/ocs/index.php/AAAI/AAAI18/paper/viewPaper/17286.
  * Hawthorne et al. (2017) Hawthorne, C., Elsen, E., Song, J., Roberts, A., Simon, I., Raffel, C., Engel, J., Oore, S., and Eck, D.  Onsets and frames: Dual-objective piano transcription.  _arXiv preprint arXiv:1710.11153_ , 2017.  URL https://arxiv.org/abs/1710.11153.
  * Holtzman et al. (2019) Holtzman, A., Buys, J., Du, L., Forbes, M., and Choi, Y.  The curious case of neural text degeneration.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=rygGQyrFvH.
  * Hsiao et al. (2021) Hsiao, W.-Y., Liu, J.-Y., Yeh, Y.-C., and Yang, Y.-H.  Compound word transformer: Learning to compose full-song music over dynamic directed hypergraphs.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pp. 178–186, 2021.  URL https://ojs.aaai.org/index.php/AAAI/article/view/16091.
  * Huang et al. (2021) Huang, J., Wang, J.-C., Smith, J. B., Song, X., and Wang, Y.  Modeling the compatibility of stem tracks to generate music mashups.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pp. 187–195, 2021.  URL https://ojs.aaai.org/index.php/AAAI/article/view/16092.
  * Jansson et al. (2017) Jansson, A., Humphrey, E., Montecchio, N., Bittner, R., Kumar, A., and Weyde, T.  Singing voice separation with deep u-net convolutional networks.  In _18th International Society for Music Information Retrieval Conference_ , pp. 23–27, 2017.  URL https://openaccess.city.ac.uk/id/eprint/19289/1/7bb8d1600fba70dd79408775cd0c37a4ff62.pdf.
  * Katharopoulos et al. (2020) Katharopoulos, A., Vyas, A., Pappas, N., and Fleuret, F.  Transformers are rnns: Fast autoregressive transformers with linear attention.  In _International Conference on Machine Learning_ , pp. 5156–5165. PMLR, 2020.  URL http://proceedings.mlr.press/v119/katharopoulos20a.html.
  * Macas et al. (2018) Macas, C., Rodrigues, A., Bernardes, G., and Machado, P.  Mixmash: A visualisation system for musical mashup creation.  In _2018 22nd International Conference Information Visualisation (IV)_ , pp. 471–477. IEEE Computer Society, 2018.  URL https://www.computer.org/csdl/proceedings-article/iv/2018/720200a471/17D45XvMcd9.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I.  Attention is all you need.  In _Advances in neural information processing systems_ , pp. 5998–6008, 2017.  URL https://proceedings.neurips.cc/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf.
  * Zheng et al. (2017) Zheng, Z., Zheng, L., and Yang, Y.  Unlabeled samples generated by gan improve the person re-identification baseline in vitro.  In _Proceedings of the IEEE international conference on computer vision_ , pp. 3754–3762, 2017.  URL https://openaccess.thecvf.com/content_iccv_2017/html/Zheng_Unlabeled_Samples_Generated_ICCV_2017_paper.html.