1

# Cache Bypassing for Machine Learning Algorithms

Asim Ikram  Department of Computer SciencesNational University of Computer and
Emerging SciencesIslamabadPakistan<EMAIL_ADDRESS>,  Muhammad Awais
Ali  Department of Computer SciencesNational University of Computer and
Emerging SciencesIslamabadPakistan<EMAIL_ADDRESS>and  Mirza Omer
Beg  Assistant ProfessorDepartment of Computer SciencesNational University of
Computer and Emerging SciencesIslamabadPakistan<EMAIL_ADDRESS>
###### Abstract.

Graphics Processing Units (GPUs) were once used solely for graphical
computation tasks but with the increase in the use of machine learning
applications, the use of GPUs to perform general purpose computing has
increased in the last few years. GPUs employ a massive amount of threads, that
in turn achieve a high amount parallelism, to perform tasks. Though GPUs have
a high amount of computation power, they face the problem of cache contention
due to the SIMT model that they use. A solution to this problem is called
”cache bypassing”. This paper presents a predictive model that analyzes the
access patterns of various machine learning algorithms and determines whether
certain data should be stored in the cache or not. It presents insights on how
well each model performs on different datasets and also shows how minimizing
the size of each model will affect its performance The performance of most of
the models were found to be around 90% with KNN performing the best but not
with the smallest size. We further increase the features by splitting the
addresses into chunks of 4 bytes. We observe that this increased the
performance of neural network substantially and increased the accuracy to
99.9% with three neurons.

cache optimization, cache bypassing, machine learning

††journalyear: 2021††doi: ††copyright: none††ccs: Software and its engineering
General programming languages††ccs: Social and professional topics History of
programming languages

##  1\. Introduction

Recent advances in artificial intelligence has led to a renewed attention
towards a diverse set of difficult combinatorial problems (Beg, 2013; Qamar et
al., 2021; Anwar and Baig, 2020; Sahar et al., 2019; Zahid et al., 2020).
Graphics Processing Units (GPUs) have been used to perform graphics intensive
tasks since the earlier days, but in the recent years, thanks to the rise of
machine learning, GPUs have been used to perform high performance tasks
(Khawaja et al., 2018). Due to the massive computational power that GPUs
offer, they are being used to perform tasks that were once meant to be
executed by CPUs. Since GPUs provide more computational power than CPUs, they
have been used extensively to train machine learning algorithms(Majeed et al.,
2020; Tariq et al., 2019; Uzair et al., 2019).

Figure 1. Miss rates for L1 cache bypassing on selected Rodinia benchmark
programs

Initially, GPUs only had a shared memory called the scratchpad memory
installed on board. While scratchpad memory was programmable and enabled rapid
fetching of data, it had its own pitfalls. The problem with scratchpad memory
was that it performed well on processes that had a uniform access pattern but
performed poorly on processes that exhibited irregular access patterns
(Huangfu and Zhang, 2015). To handle the problems that scratchpad memory
exhibited, vendors started to employ cache memory on GPUs as well.

Caches, on GPUs, perform well on tasks that exhibit non-uniform or irregular
access patterns. GPU caches do not operate in the same way as CPU caches and
hence cannot be optimized in the same way (Beg and Van Beek, 2010). They also
exhibit a low level of temporal and spatial locality. While CPU caches have
been worked upon in detail, GPU caches still remain an area to be explored (Li
et al., 2015b).

GPUs work by dispatching a large number of threads to streaming
multiprocessors (SMs) for processing large datasets in parallel. Since a large
number of threads need to be executed, the small size of the cache becomes a
bottleneck for performance especially when training deep learning models
(Naeem et al., 2020; Asad et al., 2020; Javed et al., 2020). One problem that
arises is called ”cache contention”. Under a naive approach, the processor
would evict data from the cache using the defined eviction policy when the
cache reached its maximum capacity. This approach would be ineffective since
the data already present in the cache might be more used more frequently than
newer data. Since GPUs require a massive amount of threads to be executed in
parallel on a shared cache, harmful eviction would likely occur very
frequently. One solution to this problem is called ”cache bypassing” which
involves storing some of the data directly to the L2 cache rather than the L1
cache.

Figure 1 illustrates how cache bypassing can improve performance of randomly
selected programs running on a GPU. We execute a few programs from the Rodinia
benchmark on the Kepler GPU architecture and we modify the PTX data access
instructions to randomly bypass the cache and see improvements on all cases as
seen in figure 1.

This paper proposes a mechanism that analyzes the access patterns of various
machine learning algorithms and uses a model that predicts whether a certain
address should be bypassed or not. Furthermore, since the mechanism is
intended to be embedded in the hardware, the size of the models has been
reduced to a considerable degree to reduce the cost of implementation. This
paper makes the following contributions

  * •

We propose a cache bypassing mechanism using various machine learning
algorithms that learn the access patterns of other machine learning algorithms
and takes bypassing decisions.

  * •

We shrink the size of the model and evaluate its performance. We reduce the
size of the model since it is intended to be implemented into the hardware.

  * •

We further split the addresses into multiple parts and evaluate our models on
the modified datasets and find that the performance of the neural network
increases substantially.

The rest of this paper is organized as follows. Section 2 describes the
background, section 3 presents an overview of the related work, section 4
describes our methodology, section 5 lists the experimental setup used to
conduct the experiments, section 6 lists the results of the experiments, and
section 7 summarizes our work.

##  2\. Background

###  2.1. GPU Architecture

GPUs consists of multiple parts such Streaming Multiprocessors (SMs). SMs are
processors that are designed to handle CUDA requests. Each SM has multiple
CUDA cores and a shared L1 data cache (NVIDIA Corporation, 2016). Each SM is
responsible for dispatching warps (blocks of 32 threads) to CUDA cores.

###  2.2. GPGPU

General Purpose GPU (GPGPU) are GPUs that carry out tasks that are meant to be
executed on a CPU. Since GPUs allow massive thread level parallelism, they
have been used to perform high intensive computing tasks. Due to the thread
handling capabilities of GPUs, lower end GPUs can handle tasks that might
require cutting edge CPUs.

###  2.3. Cache Bypassing

GPU caches were introduced to counteract the drawbacks of scratchpad memory.
GPU caches perform well on data that exhibits irregular access patterns but
while caches have their benefits, they also suffer from some drawbacks. GPUs
employ a mechanism called Single Instruction Multiple Threads (SIMT) that
involves multiple threads being dispatched for a single process. The small
size of the cache becomes a bottleneck when handling such a large number of
threads that run in parallel. To increase the performance of the cache in this
scenario, cache bypassing can be used. Cache bypassing involves storing only
those instructions in the cache that have a high reuse rate, the instructions
that do not have a high reuse rate a not stored in the cache and are accessed
directly from the memory. This can result in the decrease of cache miss rates.

##  3\. Related Work

GPUs have been used in the recent years for high performance computing. To
perform this tasks, vendors have started to employ caches for GPUs (Huangfu
and Zhang, 2015; Xie et al., 2015). The authors claim that due to the nature
of GPGPU applications, data data reuse rate is very low and that data can be
bypassed to improve cache performance. Yijie Huangfu and Wei Zhang (Huangfu
and Zhang, 2015) proposed a mechanism that filtered data based on the
addresses and their mechanism improved cache performance by 13.8%. The authors
in (Xie et al., 2015) classify the data into three types based on locality and
use static bypassing for data with high and low level of locality and use
dynamic bypassing for data with medium level of locality.

The authors in (Liang et al., 2015; Dai et al., 2016; Xu and Gregg, 2015)
propose bypassing mechanisms that are compiler based. The authors in (Liang et
al., 2015) state that only global load instructions are stored in the cache
and only those instructions need to be identified and bypassed (Liang et al.,
2015; Dai et al., 2016). They propose a heuristic based method for a compiler
that filters out global load instructions and generates an optimized code. The
authors in (Dai et al., 2016) propose a model that identifies the optimal
number of warps. The authors in (Xu and Gregg, 2015) propose ’Hyper Loop
Parallelism’ to improve the performance of CUDA GPUs. The authors propose a
mechanism that identifies whether a loop can be presented in a vector form or
not and build a compiler to achieve their goals. (Li et al., 2015b) also
propose a compile time framework to limiting the number of threads that can
access a cache. (Park et al., 2017) propose a bypass first policy for the last
level cache that only stores those addresses in memory that are likely to be
re-referenced. (Zhao et al., 2017) propose a bypassing scheme that is targeted
towards handling un-coalesced loads. The mechanism uses two approaches – One
is to bypass data when the number of accesses exceeds a pre-determined
threshold. The second approach is to bypass memory accesses when the L1 data
cache is stalled.

The authors in (Lee and Wu, 2016; Tian et al., 2015) propose bypassing
mechanisms that are based on reuses. (Lee and Wu, 2016) propose a mechanism
that uses feedback control loops to predict reuse patterns for each
instruction. They use a reuse table to keep track of reuses and use data from
the table to statistically determine whether to bypass an instruction or not.
The result is an almost double speedup. (Tian et al., 2015) propose an
adaptive cache bypassing mechanism to avoid premature eviction. They use the
PC trace to predict bypassing. They predict block that are likely to not be
rerefrenced and choose to bypass them which results in a higher hit rate. (Li
et al., 2015a) propose a method that dynamically bypasses instructions and
only stores those instructions in the L1 data cache that have a high reuse
rate and shore reuse distance. The authors also propose to decouple the L1
data cache to increase the energy efficiency of the cache and to enable the
storage of more reuse patterns with a lower overhead.

##  4\. Proposed Approach

The purpose of this paper is to analyze the address patterns of the machine
learning algorithms and uses a predictive model to decide whether data should
be stored in the cache or not. Since the model is intended to be embedded in
the hardware the size of the models should ideally be kept to a minimum.

A large number of machine learning models exists and using all of them was not
feasible, so we chose a subset of the algorithms available. The algorithms
that we ran our test on include Decision Tree, K Nearest Neighbors (KNN),
Logistic Regression, and Neural Network (MLP). These machine learning
algorithms consist of many parameters to tune and the tuning of these
paramaters is an exhaustive approach which is very computationally expensive.
To tune these parameters to obtain optimal values, either manually or
automatically, takes several days which was can not be considered feasible.
Therefore, we optimized these algorithms based on the core component that
builds up these algorithms. For the decision tree, we changed the depth and
impurity. We focused on these parameters because the depth controls how many
nodes exist in the tree. For increased depths, the nodes in the tree increase
along with its overall complexity. When we consider impurity, the
computational cost is higher when the impurity is lower. In the case of KNN,
we chose to change the value of K since this parameter controls the number of
clusters that the model has to construct. For logistic regression, instead of
using only one solver, we focused on trying multiple solvers to evaluate the
performance of each one and gain insights as to which one performs most
optimally on the given datasets. Neural networks have seen a massive amount of
use in the recent years. With the introduction of deep learning algorithms as
well as the increase in processing power, neural networks have been used in a
wide variety of applications such as autonomous driving and IoT (Tian et al.,
2017; Guidolini et al., 2017; Li et al., 2018; Kodali et al., 2017; Xie et
al., 2017; Mohammadi et al., 2017). Since, neural networks have been extremely
popular, we were interested in the performance of neural networks with cache
bypassing. As our model is intended to be implemented in the hardware, we
could not afford to have too many neurons in the network since that would
increase the implementation cost. We chose to change the number of neurons in
the network while keeping the number of hidden layers to one.

The dataset that we used was imbalanced and to balance it we used the SMOTE
(Chawla et al., 2002). The algorithms we used are given in algorithms 1, 2, 3,
4, 5.

for  _$DepthParam\leftarrow 1$ to $10$_ do

Initialize Decision Tree($depth=DepthParam$);

evaluate Decision Tree

end for

Algorithm 1 Decision Tree (Depth)

for  _$ImpParam\leftarrow 0$ to $0.5$_ do

Initialize Decision Tree($impurity=ImpParam$);

evaluate Decision Tree

end for

Algorithm 2 Decision Tree (Impurity)

for  _$KParam\leftarrow 1$ to $17$_ do

Initialize KNN($K=KParam$);

evaluate KNN Model

end for

Algorithm 3 KNN

for  _$SolverParam$ in $[newton-cg,LBFGS,liblinear,sag]$_ do

Initialize Logistic Regression($Solver=SolverParam$);

evaluate Logistic Regression Model

end for

Algorithm 4 Logistic Regression

for  _$NeuronParam\leftarrow 1$ to $20$_ do

Initialize NeuralNetwork($neurons=NeuronParam$);

evaluate Neural Network

end for

Algorithm 5 Neural Network

The dataset that we used consisted of only one feature that was the address.
This caused our algorithm to perform not very well, especially for the neural
network. To handle this case, we split the address into chunks of 4 bytes and
use that data to train the models once more. We present the results on for
both versions of the dataset. The algorithm that we used to split our data is
given in algorithm Splitting algorithm here

$SplitData\leftarrow[]$ ;

for  _$i$ in $range(sizeOf(Data))$_ do

$record\leftarrow data[i]$ ;

$tempList\leftarrow[]$ ;

while  _$record\not e0$_ do

$record,remainder\leftarrow record\%10$ $tempList.append(remainder)$

end while

$SplitData.append(tempList)$

end for

Algorithm 6 Splitting the Addresses

##  5\. Experimental Setup

The experiments were conducted on datasets that were generated from tensorflow
examples (Damien, 2018). The codes that we used are given in table 1. To codes
were run on the MNIST dataset (Lecun et al., 1998). We ran several iterations
on each of the examples and generated a dataset that contained approximately
1,000,000 records.

Table 1. Codes Used for Generating Datasets Code Used
---
Nearest Neighbors
Logistic Regression
Random Forest
Recurrent Neural Network

The models that we used for training are given in section 4. The models were
trained on a machine with an Intel i7 3630QM processor (2.4Ghz) 8GB RAM, and
an Nvidia Geforce GT 630M (Corporation, 2018). The specifications for the GPU
are given in table 2. The models were implemented using using Python 3.6.4. As
the model(s) would be implemented at the hardware level, we reduced the size
of the models to reduce implementation costs while maintaining an acceptable
level of accuracy.

Table 2. Geforce GT 630M Specifications Item | Value
---|---
CUDA Cores | 96
Graphics Clock | 800MHz
Memory Interface Width | 128 bits
Architecture | Kepler
Memory Bandwidth | 32 GB/s
Memory Size | 2 GB

##  6\. Results

We tested our models by using data generated from the tensorflow examples
mentioned earlier. We measured the accuracy of each model with respect to the
model’s size. The results for each of the datasets on different models are
given in the subsequent sections.

###  6.1. Case 1: One Feature

####  6.1.1. Logistic Regression Dataset

Figure 2. Results on Logistic Regression Dataset

Figure 2 shows the results obtained from the machine learning models on the
logistic regression dataset. For the decision tree, when the impurity was kept
constant and the depth of the tree was changed it was observed that the
accuracy varied when the depth was changed from 1 to 5 but from depth 5
onwards the accuracy remained the same. When the depth was kept constant and
the impurity was changed, it was observed that the accuracy was high when the
impurity was lower than 0.05 but dropped drastically when the impurity was
increased beyond this point. When the impurity was set to greater than or
equal to 0.15, the accuracy remained the same. For KNN, the accuracy decreased
more sharply in the case of uniform weights. All of the versions of the neural
networks had an accuracy of approximately 50% and exhibited an irregular
behavior when the number of neurons in the network was changed. A notable
point was that the sigmoid and lbfgs solvers gave the highest accuracy with
one neuron while the Adam solver showed a converse behavior. For logistic
regression, all the solvers gave the same error. This implies that the solver
most suited for the situation could be used without affecting the performance
of the model.

####  6.1.2. Nearest Neighbors Dataset

Figure 3. Results on the Nearest Neighbors Dataset

Figure 3 shows the results obtained on the nearest neighbors dataset using the
machine learning algorithms. In the case of the decision tree, the accuracy
increased suddenly when the depth was increased from 1 to 2 and increased very
slightly when the depth was increased from 3 to 4. The accuracy remained
constant at depths greater then or equal to 4 (impurity was constant). When
the depth was kept constant and the impurity was changed, the accuracy was
very high when the impurity was below 0.30 but decreased abruptly when the
impurity increased beyond this point. This implies that the nearest neighbors
dataset is relatively resilient to impure splitting. The accuracy obtained on
the same dataset using both versions of KNN was very high. It can be seen from
the graphs that accuracy decreased more consistently with increasing values of
K when using uniform weights. The mean absolute error obtained using logistic
regression with different solvers on the nearest neighbors dataset was the
same on all the solvers. For the neural network, the sigmoid solver exhibited
a constant accuracy when the number of neurons in the network was changed. The
LBFGS and Adam solvers displayed a completely random pattern no matter how
many neurons existed in the network. In this case, the Adam solver achieved
the highest accuracy with one neuron while LBFGS exhibited the lowest accuracy
with one neuron.

####  6.1.3. Random Forest Dataset

Figure 4. Results on the Random Forest Dataset

Figure 4 shows the results obtained on the random forest dataset using the
machine learning algorithms. I In the case of the decision tree, when the
depth was changed and impurity was kept constant, the accuracy increased at
approximately equal intervals. When the depth was kept constant and the
impurity was changed the accuracy dropped abruptly when the impurity was
increased from 0 and stayed constant at values equal to or above 0.05. This
implies that accuracy decreases substantially when the splitting of nodes is
not pure. For both of KNN, the accuracy is the highest when K is 3 but the
accuracy decreases sharply in the case of uniform weights and remains almost
constant in the case of weighted KNN. When using logistic regression, the
error on newton-cg was considerably lower then the other three solvers. The
error on LBFGS, liblinear and sag was the same. In the case of neural
networks, the sigmoid solver showed a behavior separate from the rest and
exhibited a constant accuracy no matter how many neurons existed in the
network. The LBFGS and adam solvers showed an irregular behavior when the
number of neurons was changed.

####  6.1.4. Recurrent Neural Network Dataset

Figure 5. Results on the RNN Dataset

Figure 5 shows the results obtained achieved on the RNN dataset using a the
models. For the decision tree, the accuracy increased substantially when the
depth was increased from 1 to 2 but increased by a very minor amount when the
depth was increased from 2 to 3. The accuracy remained constant at depth
greater than or equal to 3 (impurity was kept constant). When the depth was
kept constant and the impurity was increased, the accuracy decreased rapidly
when the impurity was increased from 0 meaning that a slightly impure split
can cause a rapid decrease in accuracy on this dataset. For both versions of
KNN, the highest accuracy was achieved when the value of K was 5. For values
of K greater than 5, the behavior of both versions was almost the same. The
error on all the solvers remained the same when using logistic regression.
This implies that for the RNN dataset, any of the solvers could be applied
(most preferably the best one for the situation) without affecting the error.
For the neural network, all of the solvers showed a random behavior when the
number of neurons in the network was changed. The LBFGS solver achieved the
highest accuracy when one neuron existed in the network. The sigmoid and Adam
solvers exhibited the opposite behavior under the same conditions.

###  6.2. Case 2: Multiple Features

At first we tested our model on the original dataset that consisted of only
one feature. Since using one feature was causing a linear mapping, we decided
to split the data into parts to make a more generic model and tested our
models on the modified dataset as well. The results of the modification are
described in the subsequent sections.

####  6.2.1. Logistic Regression Dataset

Figure 6. Results on the Logistic Regression Dataset (Case 2)

Figure 6 shows the results obtained from the machine learning models on the
logistic regression dataset having multiple features. For the decision tree,
when the impurity was kept constant and the depth of the tree was changed it
was observed that the accuracy varied when the depth was changed from 1 to 9
and kept increasing. When the depth was kept constant and the impurity was
changed, it was observed that the accuracy was high when the impurity was
lower than 0.05 but dropped drastically when the impurity was increased beyond
this point. When the impurity was set to greater than or equal to 0.15, the
accuracy remained the same. For KNN, the accuracy decreased more sharply in
the case of uniform weights. All of the versions of the neural networks had an
accuracy of approximately 50% and exhibited an irregular behavior when the
number of neurons in the network was changed. A notable point was that the
sigmoid and lbfgs solvers gave the highest accuracy with two neurons while the
highest accuracy was achieved with five neurons using Adam solver. For
logistic regression, all the solvers gave the same error. This implies that
the solver most suited for the situation could be used without affecting the
performance of the model.

####  6.2.2. Nearest Neighbors Dataset

Figure 7. Results on the Nearest Neighbors Dataset (Case 2)

Figure 7 shows the results obtained on the nearest neighbors dataset with
multiple features .In the case of the decision tree, the accuracy increased
suddenly when the depth was increased from 1 to 2 and increased very slightly
when the depth was increased from 3 to 4. The accuracy remained constant at
depths greater then or equal to 4 (impurity was constant). When the depth was
kept constant and the impurity was changed, the accuracy was very high when
the impurity was below 0.0.05 but decreased slightly when the impurity
increased beyond this point and then abruptly when it was increased beyond
0.25. This implies that the nearest neighbors dataset is relatively resilient
to impure splitting. The accuracy obtained on the same dataset using both
versions of KNN was very high. It can be seen from the graphs that accuracy
decreased more consistently with increasing values of K when using uniform
weights and it increased gradually with weighted KNN as the number of K was
increased. The mean absolute error obtained using logistic regression with
different solvers on the nearest neighbors dataset was the same on all the
solvers. For the neural network, the sigmoid solver exhibited a constant
accuracy when the number of neurons in the network was changed. The LBFGS and
Adam solvers displayed almost similar pattern. In LBFGS the accuracy was
initially 75% with one neuron and then there is a sudden decrease when the
number of neurons were increased to 2 but after that showed the same behavior
as Adam, it increased gradually as the number of neurons were increased.

####  6.2.3. Random Forest Dataset

Figure 8. Results on the Random Forest Dataset (Case 2)

Figure 8 shows the results obtained on the random forest dataset with multiple
features using the machine learning algorithms. In the case of the decision
tree, when the depth was changed and impurity was kept constant, the accuracy
increased gradually. When the depth was kept constant and the impurity was
changed the accuracy dropped abruptly when the impurity was increased from 0
and stayed constant at values equal to or above 0.15. This implies that
accuracy decreases substantially when the splitting of nodes is not pure. For
the KNN, the accuracy is the highest when K is 1 but the accuracy decreases
drastically in both cases when K was increased. When using logistic
regression, the error on all three solvers remain the same.In the case of
neural networks, all solvers exhibit the same behavior as the number of
neurons increased. Neural network accuracy drastically increases and then kept
constant as the neurons were further increased.

Figure 9. Results on the RNN Dataset

####  6.2.4. Recurrent Neural Network Dataset

Figure 9 shows the results achieved on the RNN dataset with multiple features.
For the decision tree, the accuracy increased substantially when the depth was
increased from 1 to 2 but then remained constant at depth greater than 2. When
the depth was kept constant and the impurity was increased, the accuracy
decreased rapidly when the impurity was increased from 0 meaning that a
slightly impure split can cause a rapid decrease in accuracy on this dataset.
For both versions of KNN, the highest accuracy was achieved when the value of
K was 5. For values of K greater than 5, the behavior of both versions was
almost the same. The error on all the solvers remained the same when using
logistic regression except for the LIBLINEAR, which shows a slight decrease in
error than the others. This implies that for the RNN dataset, LIBLINEAR
solvers could be applied (most preferably the best one for the situation)
without affecting the error. For the neural network, LBFGS shows the highest
accuracy when 3 neurons were used and then it kept constant as the number of
neurons were increased. The sigmoid and Adam solvers exhibited an irregular
behavior when the number of neurons in the network was changed.

##  7\. Conclusion

GPUs have been used in the recent years to perform general purpose computing
tasks. For this purpose, caches have been added to GPUs but they perform
poorly due to the massive number of threads accessing the cache. This paper
analyzes different machine learning algorithms and presents insights on
whether to bypass or cache addresses using various machine learning
algorithms. It also presents how reducing the size of the machine learning
algorithms will affect their performance.

## References

  * (1)
  * Anwar and Baig (2020) Talha Anwar and Omer Baig. 2020.  TAC at SemEval-2020 Task 12: Ensembling Approach for Multilingual Offensive Language Identification in Social Media. In Proceedings of the Fourteenth Workshop on Semantic Evaluation. 2177–2182.
  * Asad et al. (2020) Muhammad Asad, Muhammad Asim, Talha Javed, Mirza O Beg, Hasan Mujtaba, and Sohail Abbas. 2020\.  DeepDetect: detection of distributed denial of service attacks using deep learning.  Comput. J. 63, 7 (2020), 983–994.
  * Beg and Van Beek (2010) Mirza Beg and Peter Van Beek. 2010.  A graph theoretic approach to cache-conscious placement of data for direct mapped caches. In Proceedings of the 2010 international symposium on Memory management. 113–120.
  * Beg (2013) Mirza Omer Beg. 2013\.  Combinatorial problems in compiler optimization.  (2013).
  * Chawla et al. (2002) Nitesh V. Chawla, Kevin W. Bowyer, Lawrence O. Hall, and W. Philip Kegelmeyer. 2002.  SMOTE: Synthetic Minority Over-sampling Technique.  J. Artif. Int. Res. 16, 1 (2002), 321–357.
  * Corporation (2018) NVIDIA Corporation. 2018\.  Geforce GT 630M Specifications.  (March 2018).  Retrieved March 2, 2005 from https://www.geforce.com/hardware/notebook-gpus/geforce-gt-630m/specifications
  * Dai et al. (2016) Hongwen Dai, Chao Li, Huiyang Zhou, Saurabh Gupta, Christos Kartsaklis, and Mike Mantor. 2016\.  A Model-Driven Approach to Warp / Thread-Block Level GPU Cache Bypassing Experimental Methodology. In 53rd ACM/EDAC/IEEE Design Automation Conference (DAC). 94:1—-94:6.
  * Damien (2018) Aymeric Damien. 2018\.  TensorFlow-Examples.  (March 2018).  Retrieved March 2, 2005 from https://github.com/aymericdamien/TensorFlow-Examples
  * Guidolini et al. (2017) Ranik Guidolini, Alberto F. De Souza, Filipe Wall Mutz, and Claudine Badue. 2017.  Neural-based model predictive control for tackling steering delays of autonomous cars. In International Joint Conference on Neural Networks, IJCNN 2017,. 4324–4331.  https://doi.org/10.1109/IJCNN.2017.7966403
  * Huangfu and Zhang (2015) Yijie Huangfu and Wei Zhang. 2015.  Boosting GPU performance by profiling-based L1 data cache bypassing. In IEEE/ACM 15th International Symposium on Cluster, Cloud, and Grid Computing, CCGrid 2015. 1119–1122.  https://doi.org/10.1109/CCGrid.2015.67
  * Javed et al. (2020) Abdul Rehman Javed, Mirza Omer Beg, Muhammad Asim, Thar Baker, and Ali Hilal Al-Bayatti. 2020\.  AlphaLogger: Detecting motion-based side-channel attack using smartphone keystrokes.  Journal of Ambient Intelligence and Humanized Computing (2020), 1–14.
  * Khawaja et al. (2018) Hussain S Khawaja, Mirza O Beg, and Saira Qamar. 2018.  Domain Specific Emotion Lexicon Expansion. In 2018 14th International Conference on Emerging Technologies (ICET). 1–5.
  * Kodali et al. (2017) Sreela Kodali, Patrick Hansen, Niamh Mulholland, Paul N. Whatmough, David M. Brooks, and Gu-Yeon Wei. 2017.  Applications of Deep Neural Networks for Ultra Low Power IoT. In IEEE International Conference on Computer Design, ICCD. 589–592.  https://doi.org/10.1109/ICCD.2017.102
  * Lecun et al. (1998) Yann Lecun, Leon Bottou, Y Bengio, and Patrick Haffner. 1998\.  Gradient-Based Learning Applied to Document Recognition. In Proceedings of the IEEE, Vol. 86. 2278 – 2324.
  * Lee and Wu (2016) Shin Ying Lee and Carole Jean Wu. 2016.  Ctrl-C: Instruction-Aware Control Loop Based Adaptive Cache Bypassing for GPUs. In Proceedings of the 34th IEEE International Conference on Computer Design, ICCD 2016. 133–140.  https://doi.org/10.1109/ICCD.2016.7753271
  * Li et al. (2015b) Ang Li, Gert-Jan van den Braak, Akash Kumar, and Henk Corporaal. 2015b.  Adaptive and Transparent Cache Bypassing for GPUs. In International Conference for High Performance Computing, Networking, Storage and Analysis. 1–12.  https://doi.org/10.1145/2807591.2807606
  * Li et al. (2015a) Chao Li, Shuaiwen Leon Song, Hongwen Dai, Albert Sidelnik, Siva Kumar Sastry Hari, and Huiyang Zhou. 2015a.  Locality-Driven Dynamic GPU Cache Bypassing. In Proceedings of the 29th ACM on International Conference on Supercomputing. 67–77.  https://doi.org/10.1145/2751205.2751237
  * Li et al. (2018) He Li, Kaoru Ota, and Mianxiong Dong. 2018.  Learning IoT in Edge: Deep Learning for the Internet of Things with Edge Computing.  IEEE Network 32, 1 (2018), 96–101.
  * Liang et al. (2015) Yun Liang, Xiaolong Xie, Guangyu Sun, and Deming Chen. 2015\.  An Efficient Compiler Framework for Cache Bypassing on GPUs.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 34, 10 (2015), 1677–1690.  https://doi.org/10.1109/TCAD.2015.2424962
  * Majeed et al. (2020) Adil Majeed, Hasan Mujtaba, and Mirza Omer Beg. 2020\.  Emotion detection in Roman Urdu text using machine learning. In Proceedings of the 35th IEEE/ACM International Conference on Automated Software Engineering Workshops. 125–130.
  * Mohammadi et al. (2017) Mehdi Mohammadi, Ala I. Al-Fuqaha, Sameh Sorour, and Mohsen Guizani. 2017.  Deep Learning for IoT Big Data and Streaming Analytics: A Survey.  CoRR abs/1712.04301 (2017).
  * Naeem et al. (2020) Bilal Naeem, Aymen Khan, Mirza Omer Beg, and Hasan Mujtaba. 2020\.  A deep learning framework for clickbait detection on social area network using natural language cues.  Journal of Computational Social Science (2020), 1–13.
  * NVIDIA Corporation (2016) NVIDIA Corporation. 2016\.  NVIDIA Tesla P100 Whitepaper.  Technical Report. 45 pages.
  * Park et al. (2017) Jason Jong Kyu Park, Yongjun Park, and Scott Mahlke. 2017\.  A Bypass First Policy for Energy-Efficient Last Level Caches. In 16th International Conference on Embedded Computer Systems: Architectures, Modeling and Simulation, SAMOS. 63–70.  https://doi.org/10.1109/SAMOS.2016.7818332
  * Qamar et al. (2021) Saira Qamar, Hasan Mujtaba, Hammad Majeed, and Mirza Omer Beg. 2021.  Relationship Identification Between Conversational Agents Using Emotion Analysis.  Cognitive Computation (2021), 1–15.
  * Sahar et al. (2019) Hareem Sahar, Abdul A Bangash, and Mirza O Beg. 2019\.  Towards energy aware object-oriented development of android applications.  Sustainable Computing: Informatics and Systems 21 (2019), 28–46.
  * Tariq et al. (2019) Muhammad Tariq, Hammad Majeed, Mirza Omer Beg, Farrukh Aslam Khan, and Abdelouahid Derhab. 2019\.  Accurate detection of sitting posture activities in a secure IoT based assisted living environment.  Future Generation Computer Systems 92 (2019), 745–757.
  * Tian et al. (2017) Yuchi Tian, Kexin Pei, Suman Jana, and Baishakhi Ray. 2017\.  DeepTest: Automated Testing of Deep-Neural-Network-driven Autonomous Cars.  CoRR abs/1708.08559 (2017).
  * Tian et al. (2015) Yingying Tian, Sooraj Puthoor, Joseph L Greathouse, Bradford M Beckmann, and Daniel A Jiménez. 2015\.  Adaptive GPU cache bypassing. In 8th Workshop on General Purpose Processing using GPUs - GPGPU. 25–35.  https://doi.org/10.1145/2716282.2716283
  * Uzair et al. (2019) Ahmed Uzair, Mirza O Beg, Hasan Mujtaba, and Hammad Majeed. 2019\.  Weec: Web energy efficient computing: A machine learning approach.  Sustainable Computing: Informatics and Systems 22 (2019), 230–243.
  * Xie et al. (2015) Xiaolong Xie, Yun Liang, Yu Wang, Guangyu Sun, and Tao Wang. 2015.  Coordinated Static and Dynamic Cache Bypassing for GPUs. In IEEE 21st International Symposium on High Performance Computer Architecture, HPCA. 76–88.  https://doi.org/10.1109/HPCA.2015.7056023
  * Xie et al. (2017) Xiaofeng Xie, Di Wu, Siping Liu, and Renfa Li. 2017\.  IoT Data Analytics Using Deep Learning.  CoRR abs/1708.03854 (2017).
  * Xu and Gregg (2015) Shixiong Xu and David Gregg. 2015.  Exploiting Hyper-Loop Parallelism in Vectorization to Improve Memory Performance on CUDA GPGPU. In 14th IEEE International Conference on Trust, Security and Privacy in Computing and Communications, TrustCom. 53–60.  https://doi.org/10.1109/Trustcom.2015.612
  * Zahid et al. (2020) Rabail Zahid, Muhammad Owais Idrees, Hasan Mujtaba, and Mirza Omer Beg. 2020.  Roman Urdu reviews dataset for aspect based opinion mining. In 2020 35th IEEE/ACM International Conference on Automated Software Engineering Workshops (ASEW). 138–143.
  * Zhao et al. (2017) Chen Zhao, Fei Wang, Zhen Lin, Huiyang Zhou, and Nanning Zheng. 2017.  Selectively GPU cache bypassing for un-coalesced loads. In International Conference on Parallel and Distributed Systems - ICPADS. 908–915.  https://doi.org/10.1109/ICPADS.2016.0122


#  Mathematical structure of perfect predictive reservoir computing for
autoregressive type of time series data

Tsuyoshi Yoneda  Graduate School of Economics, Hitotsubashi University, 2-1
Naka, Kunitachi, Tokyo 186-8601, Japan<EMAIL_ADDRESS>
###### Abstract.

Reservoir Computing (RC) is a type of recursive neural network (RNN), and
there can be no doubt that the RC will be more and more widely used for
building future prediction models for time-series data, with low training
cost, high speed and high computational power. However, research into the
mathematical structure of RC neural networks has only recently begun. Bollt
(2021) clarified the necessity of the autoregressive (AR) model for gaining
the insight into the mathematical structure of RC neural networks, and
indicated that the Wold decomposition theorem is the milestone for
understanding of these. Keeping this celebrated result in mind, in this paper,
we clarify hidden structures of input and recurrent weight matrices in RC
neural networks, and show that such structures attain perfect prediction for
the AR type of time series data.

###### Key words and phrases:

Reservoir computing, Autoregressive model, universal approximation theorem,
almost periodic functions, transcendental numbers

###### 2020 Mathematics Subject Classification:

Primary 68T27; Secondary 11B50, Tertiary 42A16

##  1\.  Introduction

Reservoir Computing (RC) is a type of recursive neural network (RNN). Gilpin
[4] evaluated 24 statistical forecasting models across 135 dynamical systems,
including RC, autoregressive moving averages (ARIMA), deep neural networks
such as the transformer model, long-short-term-memory networks (LSTM), vanilla
recurrent neural networks (RNN), temporal convolutional neural networks and
neural basis expansion/neural hierarchical interpolation (NBEATS/NHiTS). The
best-performing machine learning models require very long training times, in
contrast, the RC exhibits competitive performance with two orders of magnitude
less training time. Thus there can be no doubt that the RC will be more and
more widely used for building future prediction models for time-series data,
with low training cost, high speed and high computational power.

On the other hand, research into the mathematical structure of RC neural
networks has only recently begun. Bollt [1] clarified the necessity of the
autoregressive (AR) model for gaining the insight into the mathematical
structure of RC neural networks, and indicated the Wold decomposition theorem
[10] is the milestone for understanding of these. More precisely, in the
stochastic framework, a zero mean covariance stationary vector process admits
a vector AR representation (see [1, Section V]). Furthermore, Gauthier et.al
[3] proposed a next generation RC with quadratic reservoir vectors, which
focuses not only on the mathematical understanding of the RC, but also on the
fundamental improvement of it. In contrast to these celebrated results, we
stick to the deterministic framework, and in this paper, we clarify hidden
structures of input and recurrent weight matrices, and show that these
structures attain perfect prediction for the AR type of time series data.

##  2\. AR model and almost periodic functions

Before going into structures of input and recurrent weight matrices, first we
construct both training and reference data. Let us start from the following
condition on smooth functions $\phi\in C^{\infty}(\mathbb{R})\cap
L^{\infty}(\mathbb{R})$, which is naturally expressing “recurring pattern”:

(1) $\begin{cases}\text{From any sequence of the form
$\\{\phi(t+h_{n})\\}_{n}$ where $h_{n}$ are real numbers},\\\ \text{one can
extract a subsequence converging uniformly on the real line.}\end{cases}$

Due to Corduneanu [2, Theorems 1.10 and 1.11], “recurring pattern” is nothing
more than the almost periodicity (necessary and sufficient conditions),
expressed as follows:

(2)
$\phi(t)=\sum_{\lambda\in\Lambda}a_{\lambda}\sin\left(\lambda(t-b_{\lambda})\right),\quad\\{a_{\lambda}\\}_{\lambda},\\{b_{\lambda}\\}_{\lambda}\subset\mathbb{R},\quad\Lambda(\subset\mathbb{R})\text{
is countable.}$

######  Remark 1.

We see that almost periodic functions possess quasi-periodic orbits, so, these
are integrable systems (see the well-known Arnold-Liouville theorem). We now
explain it briefly. Let $L\in\mathbb{Z}_{\geq 1}$,
$\\{\lambda_{j}\\}_{j=1}^{L}\subset\mathbb{R}$ and let $\mathcal{M}$ be a
torus such that

$\mathcal{M}=\prod_{j=1}^{L}(\mathbb{R}/(2\pi\mathbb{Z})).$

Also let $x_{t}$ be a shift operator (diffeomorphism) such that

$x_{t}=x_{0}+\tau t:\mathcal{M}\to\mathcal{M}\quad(x_{0}\mapsto x_{t}),\quad
t\in\mathbb{R},\quad\tau=\\{\lambda_{j}\\}_{j=1}^{L}\in
T_{x_{t}}\mathcal{M}\cong\mathbb{R}^{L}.$

Then there exists a $g:\mathcal{M}\to\mathbb{R}$ such that

$\phi(t)=\sum_{j=1}^{L}a_{\lambda_{j}}\sin\left(\lambda_{j}(t-b_{\lambda_{j}})\right)=g\circ
x_{t}.$

More specifically, we set $g$ as follows:

$g(t_{1},t_{2},\cdots,t_{L})=\sum_{j=1}^{L}a_{\lambda_{j}}\sin t_{j}.$

This expression exhibits nothing more than quasi-periodic orbit. Kobayashi et.
al. [8] investigated the RC from a dynamical system perspective, such as
unstable fixed points, periodic orbits, chaotic saddle, Lyapunov exponents and
manifold structures (see also [6, 9]). We see that their unstable periodic
orbit must be related to our quasi-periodic orbit, since the definition of
“chaos” needs a sort of the following property (see [7]):

  * •

Let $f$ be a map, which takes an interval $I$ to itself. Then periodic points
of $f$ are dense in $I$.

We emphasize that almost periodic functions are indispensable for
mathematically analyzing the AR model (in the deterministic framework). For
$\\{p_{\ell}\\}_{\ell=1}^{L}\subset\mathbb{R}$, the AR model is described as
follows:

(3)
$y(t)=p_{1}y(t-1)+p_{2}y(t-2)+\cdots+p_{L-1}y(t-L+1)+p_{L}y(t-L)\quad(t\geq
0)$

with prescribed initial data $\\{y(-\ell)\\}_{\ell=1}^{L}$. We now explain
that this AR model crucially includes the structure of almost periodic
functions (2). We plug the following initial data (looking into the
characteristic equation)

$y(-\ell)=\mu^{L-\ell},\quad(\mu\in\mathbb{R},\ \ell=0,1,\cdots,L)$

into (3). Throughout this paper, we choose $L>0$ to be even, and we look into
eigenfunctions of the characteristic equation whose modulus of eigenvalues are
exactly $1$. In this context, the following equality is the crucial:

$\begin{split}(\mu-e^{i\lambda})(\mu-e^{-i\lambda})&=\mu^{2}-2\mu\cos\lambda+1\quad\text{for}\quad\mu,\lambda\in\mathbb{R}.\end{split}$

We multiply this type of second order polynomials $L/2$ times, then we obtain
the following equality which clarifies the relation between almost periodicity
and the AR model (factorization of $L$-th degree polynomial):

$0=-\sum_{\ell=1}^{L}p_{\ell}\mu^{L-\ell}+\mu^{L}=\prod_{j=1}^{L/2}(\mu^{2}-2\cos\lambda_{j}\mu+1),$

namely, if $\\{p_{\ell}\\}_{\ell}$ satisfies the above equality (at least we
can easily figure out that $p_{L}=-1$), $\\{y(t)\\}_{t\geq 0}$ for the AR
model (3) can be expressed as follows:

$y(t)=\sum_{j=1}^{L/2}a_{j}\sin\left(\lambda_{j}(t-b_{j})\right),\quad
t=0,1,\cdots,$

where $a_{j},b_{j}\in\mathbb{R}$ are uniquely determined by the initial data.
Since the almost periodic functions naturally possess recurring pattern, in
the next section, we employ this AR model data (almost periodic functions) as
the both training and reference data, more precisely,

  * •

$\\{y(t)\\}_{t\in\mathbb{Z}_{<-L}}$ as the training data,

  * •

$\\{y(t)\\}_{t\in\mathbb{Z}_{\geq-L}}$ as the reference data.

##  3\. Mathematical structure of reservoir computing (main result)

Throughout this paper we set $N$ as the number of reservoir nodes (we will
determine this $N$ later). First we formulate the RC and then we state the
main theorem. For sufficiently small $\varepsilon>0$, let
$h:\mathbb{R}\to[-1,1]$ be an activate function (which is allowed to be odd
symmetric) as follows:

(4) $|h(t)-\tanh t|<\varepsilon\quad\text{for}\quad t\in\mathbb{R}.$

However, we expect that this condition (4) is not needed, that is, the setting
$h(t)=\tanh t$ may goes well. To give a simpler proof of the mathematical main
theorem, we decided to employ (4). Now let us discretize the range $[-1,1]$ as
follows: For $K\in\mathbb{Z}_{\geq 1}$, we choose $\\{a_{k}^{K}\\}_{k=0}^{2K}$
such that (we employ transcendental numbers)

  * •

$\\{a_{1}^{K},a_{2}^{K},\cdots,a^{K}_{K-1},a_{K+1}^{K},\cdots,a_{2K-1}^{K}\\}\subset\\{\pm
e^{-\frac{n}{m}}\ ;m,n\in\mathbb{Z}_{\geq
1}\\}\subset\mathbb{R}\setminus\mathbb{Q}$,

  * •

$-1=a_{0}^{K}<a_{1}^{K}<a_{2}^{K}<\cdots<a_{K}^{K}=0<a_{K+1}^{K}<\cdots<a_{2K-1}^{K}<a_{2K}^{K}=1$,

  * •

$\displaystyle\lim_{K\to\infty}\sup_{1\leq k\leq
2K}|a_{k-1}^{K}-a_{k}^{K}|=0$.

By the Lindemann–Weierstrass theorem, we see that

(5)
$\displaystyle\frac{a_{k^{\prime}}^{K}}{a_{k}^{K}}\in(\mathbb{R}\setminus\mathbb{Q})\cup\\{0\\}\cup\\{-1\\}\quad(k\not=k^{\prime},\
k\not=K,\ k,k^{\prime}\geq 1),$

(6)
$\displaystyle\sum_{\ell=1}^{L}\frac{a_{k^{\prime}_{\ell}}^{K}}{a_{k_{\ell}}^{K}}\in(\mathbb{R}\setminus\mathbb{Q})\cup\\{-L,-L+1,\cdots,-1,0,1,\cdots,L-1\\}\quad(k_{\ell}\not=K,\
k_{\ell},k^{\prime}_{\ell}\geq 1),$

except for the $k_{1}=k^{\prime}_{1},\
k_{2}=k^{\prime}_{2},\cdots,k_{L}=k^{\prime}_{L}$ case.

######  Remark 2.

$\displaystyle\sum_{\ell=1}^{L}\frac{a_{k^{\prime}_{\ell}}^{K}}{a_{k_{\ell}}^{K}}=L$

if and only if $k_{1}=k^{\prime}_{1},\
k_{2}=k^{\prime}_{2},\cdots,k_{L}=k^{\prime}_{L}$.

In what follows, we employ the AR model data (almost periodic functions) as
the both training and reference data:

(7)
$y(t)=\sum_{\ell=1}^{L}a_{\ell}\sin(\lambda_{\ell}(t-b_{\ell}))=\sum_{\ell=1}^{L}p_{\ell}y(t-\ell),\quad
t\in\mathbb{Z},$

for some suitably prescribed $\\{p_{\ell}\\}_{\ell=1}^{L}$,
$\\{\lambda_{\ell}\\}_{\ell=1}^{L}$, $\\{a_{\ell}\\}_{\ell=1}^{L}$ and
$\\{b_{\ell}\\}_{\ell=1}^{L}$, with the normalization $y(t)\in[-1,1]$
($t\in\mathbb{Z}$). We now discretize this $y(t)$. There exists a unique
$k_{t}\in\\{1,\cdots,2K\\}$ such that

(8)
$\begin{split}\frac{a_{k_{t}-1}^{K}+a_{k_{t}}^{K}}{2}<&y(t)\leq\frac{a_{k_{t}}^{K}+a^{K}_{k_{t}+1}}{2}\quad(k_{t}=2,3,\cdots,2K-2)\quad\text{or}\\\
a_{0}^{K}<&y(t)\leq\frac{a_{1}^{K}+a^{K}_{2}}{2}\quad(k_{t}=1)\quad\text{or}\\\
\frac{a_{2K-2}^{K}+a_{2K-1}^{K}}{2}<&y(t)\leq a_{2K}^{K}\quad(k_{t}=2K-1),\\\
\end{split}$

thus we can appropriately define the discretized $\bar{y}$ as follows:

(9) $\bar{y}(t):=a_{k_{t}}^{K}.$

Note that, we can simplify this discretization as follows:

$\bar{y}(t)=\operatorname*{arg\,min}_{a\in\\{a^{K}_{k}\\}_{k=1}^{2K-1}}|y(t)-(a-0)|,$

where $a-0:=a-\varepsilon$ for any sufficiently small $\varepsilon>0$. Now we
determine the training time steps $T>0$. To determine it, we just apply (1),
namely, there exists a sufficiently large $T>0$ such that

(10) $\sup_{t}|y(t-T)-y(t)|\ll 1/K.$

This means that the sequence pattern

$\bar{y}(-L),\bar{y}(-L+1)\cdots,\bar{y}(-1)$

is almost the same as

$\bar{y}(-L-T),\bar{y}(-L+1-T),\cdots,\bar{y}(-1-T).$

Rigorously, it may still have an error $\ll 1/K$, but for simplicity, here, we
identify these two sequences. Now we set up the RC as follows:

  * •

Training phase

From time-series (training) data
$\\{\bar{y}(t)\\}_{t\in\mathbb{Z}_{<-L}\cap\mathbb{Z}_{\geq-T}}$, we create
reservoir state vectors (column vectors)
$\\{\bar{r}(t)\\}_{t\in\mathbb{Z}_{<-L}\cap\mathbb{Z}_{\geq-T}}\subset\mathbb{R}^{N}$
by using the following RC: For each fixed $t$, first, we determine the
following tentative reservoir state vectors $\widetilde{r}_{t}(t-\ell+1)$
($\ell=L,L-1,\cdots,2,1$ in this order) inductively:

(11) $\begin{split}\widetilde{r}_{t}(t-L+1)&=h(W^{in}\bar{y}(t-L)),\\\
\widetilde{r}_{t}(t-\ell+1)&=h(W\widetilde{r}_{t}(t-\ell)+W^{in}\bar{y}(t-\ell))\quad(\ell=L-1,L-2,\cdots,2,1),\end{split}$

where $W^{in}\in\mathbb{R}^{N\times 1}$ is a column vector (degenerated input
weight matrix), $W\in\mathbb{R}^{N\times N}$ is a square matrix (recurrent
weight matrix). These $W^{in}$ and $W$ are prescribed vector and matrix, and
we will explain concrete $W^{in}$ and $W$ in the next section. Then we set

$\bar{r}(t):=\widetilde{r}_{t}(t).$

Note that, this $L$ should be corresponding to the transient time interval in
the usual RC.

######  Remark 3.

In this paper, we neglect the term $\widetilde{r}_{t}(t-L)$ (in the first
equation in (11)) which exists in the usual RC. Even if we take
$\widetilde{r}_{t}(t-L)$ into account, the contribution of this term is
relatively small, if the recurrent weight matrix $W$ satisfies the echo state
property (see Jaeger [5]), however, it remains open question whether or not
$W$ in Theorem 1 really satisfies it, when $K$ and $L$ are relatively large.

From reservoir state vectors, we determine a row vector
$W^{out}\in\mathbb{R}^{1\times N}$ (degenerated output weight matrix) by using
the mean-square error. More precisely, we find $W^{out}$ such that

(12)
$W^{out}:=\operatorname*{arg\,min}_{\widetilde{W}^{out}}\sum_{t=-T}^{-L-1}\left|y(t)-\widetilde{W}^{out}\bar{r}(t)\right|^{2}.$

  * •

Inference phase

We plug $W^{in}$, $W$ and $W^{out}$ into the following RC, and create a series
of future prediction $\\{\bar{u}(t)\\}_{t\geq 0}$ from initial reference data
$\\{\bar{y}(-\ell)\\}_{\ell=1}^{L}$:

(13) $\begin{split}&\begin{cases}\bar{r}(-L+1)&=h(W^{in}\bar{y}(-L)),\\\
\bar{r}(-\ell+1)&=h(W\bar{r}(-\ell)+W^{in}\bar{y}(-\ell)),\\\
\bar{u}(0)&=W^{out}\bar{r}(0)-\bar{\delta}_{n_{0}},\end{cases}\quad(\ell=L-1,L-2,\cdots
2,1)\\\
&\begin{cases}\bar{r}(t)&=h(W\bar{r}(t-1)+W^{in}\bar{u}(t-1)-W^{in}\bar{u}(t-L-1)),\\\
\bar{u}(t)&=W^{out}\bar{r}(t)-\bar{\delta}_{n_{t}}.\end{cases}\quad(t=1,2,\cdots).\end{split}$

where $\bar{\delta}_{n}$ is defined in (16) as averages of the errors
$\\{\delta(t)\\}_{t\in\mathbb{Z}_{<-L}\cap\mathbb{Z}_{\geq-T}}$ in Remark 6,
and the index $n_{t}$ ($t=0,1,\cdots$) is uniquely determined. See Remark 7.

######  Remark 4.

Since we do not know whether or not this $W$ possesses the echo state property
so far, thus we need to subtract $W^{in}\bar{u}(t-L-1)$ and eliminate the past
contribution.

Then we can state the main theorem as follows:

######  Theorem 1.

(Perfect prediction for $K\to\infty$.) For each $K\in\mathbb{Z}_{\geq 1}$,
there exist $h$ with (4), $W$ and $W^{in}$ such that

(14) $|\bar{u}(t)-y(t)|\lesssim_{L}\frac{2^{t}}{K}\quad(t\geq 0).$

######  Remark 5.

In the Fourier analysis, existence of the perfect prediction is rather obvious
due to (1). The point is that we found an explicit representation (i.e.
pattern memory) of it.

##  4\. Proof of main theorem

The crucial point of the proof is to construct suitable $W$ and $W^{in}$. In
order to do so, we need to define a row vector which represents
$N$-consecutive time series data:

$V_{\ell}:=(V_{\ell,1},V_{\ell,2},\cdots,V_{\ell,N})\quad\text{for}\quad\ell=1,2,\cdots,L.$

Let $a_{k}:=a_{k}^{K}$. First let $\sigma_{j}$ ($j=1,2,\cdots,N$) be a
permutation operator, namely,

$\sigma_{j}:\\{1,2,\cdots,L\\}\to\\{a_{1},a_{2},\cdots,a_{K-1},0,a_{K+1},a_{K+2},\cdots,a_{2K-1}\\}$

($\ell\mapsto\sigma_{j}(\ell)$) and $\sigma_{j}\not=\sigma_{j^{\prime}}$
($j\not=j^{\prime}$). We exclude the case when

$\sigma_{j}(\ell)=0\quad\text{for}\quad\ell=1,2,\cdots,L\quad(\text{identically
zero case}),$

and we impose the following two conditions for uniquely determining $N$:

  * •

For any $t\in\mathbb{Z}_{<-L}\cap\mathbb{Z}_{\geq-T}$, there is
$j\in\\{1,\cdots,N\\}$ such that $\sigma_{j}(\ell)=\bar{y}(t-\ell)$ for
$\ell=1,2,\cdots,L$,

  * •

For any $j\in\\{1,\cdots,N\\}$ there is
$t\in\mathbb{Z}_{<-L}\cap\mathbb{Z}_{\geq-T}$ such that
$\sigma_{j}(\ell)=\bar{y}(t-\ell)$ for $\ell=1,2,\cdots,L$.

Note that $N\leq(2K-1)^{L}-1$ due to the sequence with repetition. Then we can
define the representation of $L$-consecutive time series data as follows:

$V_{\ell,j}:=\sigma_{j}(\ell).$

This definition covers all patterns of $L$-consecutive time series data in
$\\{\bar{y}(t)\\}_{t\in\mathbb{Z}_{<-L}\cap\mathbb{Z}_{\geq-T}}$, in other
words, in the training phase, there exists a column vector

$e:=(\underbrace{0,0,\cdots,0,1,0\cdots,0}_{N})^{T}$

such that

$V_{\ell}e=\bar{y}(t-\ell)\quad(\ell=1,\cdots,L).$

In particular, due to (10), even for the initial reference data
$\\{\bar{y}(-\ell)\\}_{\ell=1}^{L}$, there exists a column vector

$e:=(\underbrace{0,0,\cdots,0,1,0\cdots,0}_{N})^{T}$

such that

$V_{\ell}e=\bar{y}(-\ell)\quad(\ell=1,\cdots,L).$

In order to obtain the future prediction data
$\\{\bar{u}(t)\\}_{t=0}^{\infty}$ sequentially, we need to classify each
patterns as follows:

(15)
$\mathcal{T}_{n}:=\left\\{t\in\mathbb{Z}_{<-L}\cap\mathbb{Z}_{\geq-T}:\sigma_{n}(\ell)=\bar{y}(t-\ell)\quad\text{for}\quad\ell=1,2,\cdots,L\right\\}.$

The crucial point is that, for $t\in\mathcal{T}_{n}$,

$\operatorname*{arg\,min}_{a\in\\{a_{k}^{K}\\}_{k=1}^{2K-1}}\left|\sum_{\ell=1}^{L}p_{\ell}y(t-\ell)-(a-0)\right|$

may NOT be uniquely determined. In this case we just choose arbitrary one
$t^{*}\in\mathcal{T}_{n}$ such that

$v^{*}_{n}:=\operatorname*{arg\,min}_{a\in\\{a_{k}^{K}\\}_{k=1}^{2K-1}}\left|\sum_{\ell=1}^{L}p_{\ell}y(t^{*}-\ell)-(a-0)\right|$

and we define the row vector $V^{*}\in\mathbb{R}^{N}$ as follows:

$V^{*}:=(v^{*}_{1},v^{*}_{2},\cdots,v^{*}_{N}).$

We need this $V^{*}$ in the inference phase, to obtain the future prediction
data $\\{\bar{u}(t)\\}_{t=0}^{\infty}$ sequentially.

######  Remark 6.

We observe the controllable errors:

$\begin{split}&v^{*}_{n}-\sum_{\ell=1}^{L}p_{\ell}\bar{y}(t-\ell)=:\delta^{*}(t)\quad\text{for}\quad
t\in\mathcal{T}_{n},\\\
&|\delta^{*}(t)|\lesssim\left|\sum_{\ell=1}^{L}p_{\ell}y(t^{*}-\ell)-\sum_{\ell=1}^{L}p_{\ell}\bar{y}(t-\ell)\right|+\frac{1}{K}\lesssim_{L}\frac{1}{K},\\\
&y(t)-\sum_{\ell=1}^{L}p_{\ell}\bar{y}(t-\ell)=\sum_{\ell=1}^{L}p_{\ell}y(t-\ell)-\sum_{\ell=1}^{L}p_{\ell}\bar{y}(t-\ell)=:\delta(t),\\\
&|\delta(t)|\lesssim_{L}\frac{1}{K}.\end{split}$

Let $\bar{\delta}=(\bar{\delta}_{1},\bar{\delta}_{2},\cdots,\bar{\delta}_{N})$
be the corresponding mean averages:

(16)
$\begin{split}\bar{\delta}_{n}:=&\frac{1}{|\mathcal{T}_{n}|}\sum_{t\in\mathcal{T}_{n}}(\delta(t)-\delta^{*}(t)),\quad\text{in
other words,}\\\
\bar{\delta}_{n}:=&\operatorname*{arg\,min}_{\delta}\frac{1}{|\mathcal{T}_{n}|}\sum_{t\in\mathcal{T}_{n}}|\delta(t)-\delta^{*}(t)-\delta|^{2}.\end{split}$

Clearly, $|\bar{\delta}_{n}|\lesssim_{L}1/K$.

######  Remark 7.

In the inference phase, for each $t\geq 0$, there exists a unique column
vector

$e_{n_{t}}:=(\underbrace{\overbrace{0,0,\cdots,0,1}^{n_{t}},0\cdots,0}_{N})^{T}$

such that

(17) $\bar{u}(t-\ell)=V_{\ell}e_{n_{t}}\quad\text{for}\quad\ell=1,2,\cdots,L.$

We have used this index $n_{t}$ in (13).

Next we define column vectors (conjugate type of vectors)

$W^{in}_{\ell}:=(W^{in}_{\ell,1},W^{in}_{\ell,2},\cdots,W^{in}_{\ell,N})^{T}\quad\text{for}\quad\ell=1,2,\cdots,L,$

as follows: First let $\sigma^{*}_{i}$ ($i=1,2,\cdots,N$) be an adjoint type
of permutation operator, namely, let

$\begin{cases}\displaystyle\sigma^{*}_{j}(\ell):=\frac{1}{\sigma_{j}(\ell)}&\quad\text{if}\quad\sigma_{j}(\ell)\not=0,\\\
\sigma^{*}_{j}(\ell):=0&\quad\text{if}\quad\sigma_{j}(\ell)=0\end{cases}$

for $\ell\in\\{1,\cdots,L\\}$ and $j\in\\{1,\cdots N\\}$. Then we can define
the conjugate type of representation of $L$-consecutive time series as
follows:

(18)
$W^{in}_{\ell,i}:=\sigma^{*}_{i}(\ell)\times\frac{1}{\\#\\{\sigma^{*}_{i}(\ell)\not=0;\ell=1,2,\cdots,L\\}}.$

For notational convenience, we set $W^{in}_{0}:=W^{in}_{L}$, also, let
$W^{in}:=W^{in}_{L}$. By the definition of this $\\{W^{in}_{\ell}\\}_{\ell}$,
then we can construct a suitable matrix $W$. More precisely, our main task now
is to assure existence of the inverse of matrix $X$:

$X:=h(W^{in}_{0}V_{1}+W^{in}_{1}V_{2}+\cdots+W^{in}_{L-1}V_{L}).$

Note that, by using this expression, the RC (11) can be rewritten as follows:

(19) $WX=W^{in}_{1}V_{1}+W^{in}_{2}V_{2}+\cdots+W^{in}_{L}V_{L}=:Y.$

######  Lemma 2.

$X$ is a regular matrix. In other words, we have $W=YX^{-1}$.

######  Remark 8.

By using this $W$, we have

$\begin{split}Wh(W^{in}\bar{y}(t-L))&=W^{in}_{1}\bar{y}(t-L),\\\
W(h(W^{in}\bar{y}(t-L+1)+W^{in}_{1}\bar{y}(t-L)))&=W^{in}_{1}\bar{y}(t-L+1)+W^{in}_{2}\bar{y}(t-L),\\\
&\vdots\\\
W\left(h\left(\sum_{\ell=1}^{L}W^{in}_{\ell-1}\bar{y}(t-\ell)\right)\right)&=\sum_{\ell=1}^{L}W^{in}_{\ell}\bar{y}(t-\ell).\\\
\end{split}$

Note that, if $\bar{y}(t-L)=0$, then we just skip the first step and start
from the second step. If $\bar{y}(t-L)=\bar{y}(t-L+1)=0$, then we skip the
first and second steps and start from the third step, and so on.

###### Proof.

The key ingredient of the proof is the following:

  * •

Let $f$ be a non-zero polynomial in $N$ variables. Then the complement of the
zero point set, that is, $\\{x\in\mathbb{R}^{N};f(x)\not=0\\}$ is dense in
$\mathbb{R}^{N}$.

In the following proof, we give a concrete representation of such density. By
(5) and (18) (see also Remark 2), we see

$\begin{split}W^{in}_{\ell-1,i}V_{\ell,j}&\in\bigg{\\{}1,2,\cdots,\\#\\{\sigma_{i}(\ell);\ell=1,2,\cdots,L\\}\bigg{\\}}\times\frac{1}{\\#\\{\sigma_{i}(\ell);\ell=1,2,\cdots,L\\}}\\\
\quad\text{or}\quad&\in(\mathbb{R}\setminus\mathbb{Q})\cup\\{0\\}.\end{split}$

By (6), Remark 2 and Lindemann–Weierstrass theorem, we see

$\sum_{\ell=1}^{L}W^{in}_{\ell-1,j}V_{\ell,j}=1\quad\text{and}\quad\sum_{\ell=1}^{L}W^{in}_{\ell-1,i}V_{\ell,j}\not=1$

for $i\not=j$. In order to construct an appropriate $h$, we use a finite set
$G$ as follows:

$G:=\left\\{\sum_{\ell=1}^{L}W^{in}_{\ell-1,i}V_{\ell,j};i,j\in\\{1,2,\cdots,N\\}\right\\}\subset\mathbb{R}.$

Note that $1\in G$. Now we take a smooth function $h:\mathbb{R}\to[-1,1]$
satisfying the following:

$\begin{split}&h(1)\in\mathbb{Q}\setminus\\{0\\},\\\ &h(\gamma)\in\left(\\{\pm
e^{-\frac{n}{m}}\ ;m,n\in\mathbb{Z}_{\geq
1}\\}\cup\\{0\\}\right)\subset(\mathbb{R}\setminus\mathbb{Q})\cup\\{0\\},\quad\text{for}\quad\gamma\subset
G\setminus\\{1\\}.\end{split}$

Then we can easily check that (applying the Lindemann–Weierstrass theorem)

  * •

$h(\gamma_{1})h(\gamma_{2})\cdots
h(\gamma_{N})\in(\mathbb{R}\setminus\mathbb{Q})\cup\\{0\\}$ for
$\\{\gamma_{n}\\}_{n=1}^{N}\in
G^{N}\setminus\\{\underbrace{1,1,\cdots,1}_{N}\\}$,

  * •

for any $\\{\tau_{n^{\prime}}\\}_{n^{\prime}=1}^{N!-1}\in\\{-1,1\\}^{N!-1}$
and $\\{\\{\gamma_{n,n^{\prime}}\\}_{n=1}^{N}\\}_{n^{\prime}=1}^{N!-1}\subset
G^{N}\setminus\\{\underbrace{1,1,\cdots,1}_{N}\\}$,

$\displaystyle\sum_{n^{\prime}=1}^{N!-1}\tau_{n^{\prime}}h(\gamma_{1,n^{\prime}})h(\gamma_{2,n^{\prime}})\cdots
h(\gamma_{N,n^{\prime}})\in(\mathbb{R}\setminus\mathbb{Q})\cup\\{0\\}.$

By applying the above properties, then we see that the determinant of the
matrix $X$ is nonzero, since it is expressed as

$|X|=\eta_{1}+\eta_{2},\quad\eta_{1}\in\mathbb{Q}\setminus\\{0\\},\quad\eta_{2}\in(\mathbb{R}\setminus\mathbb{Q})\cup\\{0\\}.$

∎

Now we resume the proof of the main theorem. First we solve the following:

$\bar{\delta}+V^{*}=W^{out}X$

Since there exists the inverse matrix $X^{-1}$, we obtin

(20) $\left(\bar{\delta}+V^{*}\right)X^{-1}=W^{out}.$

We now check that this is the desired $W^{out}$. By Remark 6 and (16), we have
that

$\begin{split}&\sum_{t}\left|y(t)-W^{out}\bar{r}(t)\right|^{2}\\\
=&\sum_{t}\left|\delta(t)-\delta^{*}(t)+V^{*}-W^{out}\bar{r}(t)\right|^{2}\\\
=&\sum_{n=1}^{N}\sum_{t\in\mathcal{T}_{n}}|\delta(t)-\delta^{*}(t)-\bar{\delta}_{n}|^{2}\\\
&\quad+2\sum_{n=1}^{N}\sum_{t\in\mathcal{T}_{n}}(\delta(t)-\delta^{*}(t)-\bar{\delta}_{n})\left(\bar{\delta}_{n}+V^{*}e_{n}-W^{out}h(Y)e_{n}\right)\\\
&\quad+\sum_{n=1}^{N}\sum_{t\in\mathcal{T}_{n}}\left|\bar{\delta}_{n}+V^{*}e_{n}-W^{out}h(Y)e_{n}\right|^{2},\end{split}$

where $e_{n}$ is a suitable column vector such that

$e_{n}:=(\underbrace{\overbrace{0,0,\cdots,0,1}^{n},0\cdots,0}_{N})^{T}.$

Therefore the minimum value of (12) can be attained by (16) and (20), which is
zero. In the inference phase (13), we show (14). First we estimate the case
$t=0$. Since

$\bar{u}(0)=v^{*}_{n_{0}}=\sum_{\ell=1}^{L}p_{\ell}\bar{y}(-\ell)+\delta^{*}(0),\\\
\quad\text{and}\quad y(0)=\sum_{\ell=1}^{L}p_{\ell}\bar{y}(-\ell)+\delta(0),$

we have

$|\bar{u}(0)-y(0)|\lesssim|\delta(0)|+|\delta^{*}(0)|\lesssim_{L}\frac{2}{K}.$

Next we estimate the case $t=1$. Since

$\begin{split}\bar{u}(1)=v^{*}_{n_{1}}&=p_{1}\bar{u}(0)+\sum_{\ell=2}^{L}p_{\ell}\bar{y}(-\ell+1)+\delta^{*}(1)\\\
&=p_{1}(y(0)+\delta^{*}(0)-\delta(0))+\sum_{\ell=2}^{L}p_{\ell}\bar{y}(-\ell+1)+\delta^{*}(1)\\\
\end{split}$

and

$\begin{split}y(1)&=\sum_{\ell=1}^{L}y(1-\ell)=\sum_{\ell=1}^{L}\bar{y}(1-\ell)+\delta(1),\end{split}$

we have

$\begin{split}|\bar{u}(1)-y(1)|\lesssim_{L}\frac{4}{K}.\end{split}$

Also we estimate the case $t=2$. Since

$\begin{split}\bar{u}(2)=v^{*}_{n_{2}}&=p_{1}\bar{u}(1)+p_{2}\bar{u}(0)+\sum_{\ell=3}^{L}p_{\ell}\bar{y}(-\ell+2)+\delta^{*}(2)\\\
&=p_{1}(y(1)+\delta^{*}(1)-\delta(1)+\delta^{*}(0)-\delta(0))+p_{2}(y(0)+\delta^{*}(0)-\delta(0))\\\
&\quad+\sum_{\ell=3}^{L}p_{\ell}\bar{y}(-\ell+2)+\delta^{*}(2)\\\ \end{split}$

and

$\begin{split}y(2)&=\sum_{\ell=1}^{L}y(2-\ell)=\sum_{\ell=1}^{L}\bar{y}(2-\ell)+\delta(2),\end{split}$

we have

$\begin{split}|\bar{u}(2)-y(2)|\lesssim_{L}\frac{8}{K}.\end{split}$

Repeating this argument, we have

$|\bar{u}(t)-y(t)|\lesssim_{L}\frac{2^{t}}{K}.$

This is the desired estimate.

Acknowledgments. The author is grateful to Professors Chikara Nakayama and
Yoshitaka Saiki for valuable comments. Research of TY was partly supported by
the JSPS Grants-in-Aid for Scientific Research 20H01819.

Conflict of Interest. The authors have no conflicts to disclose.

## References

  * [1] Bollt, E. (2021). On explaining the surprising success of reservoir computing forecaster of chaos? The universal machine learning dynamical system with contrast to VAR and DMD. Chaos, 31, 013108.
  * [2] Corduneanu, C. (1968). Almost Periodic Functions. Interscience Publishers, New York.
  * [3] Gauthier, D. J., Bollt, E., Griffith. A. & Barbosa, W. A. S. (2021). Next generation reservoir computing. Nature Communications, 12, 5564.
  * [4] W. Gilpin, W. (2023). Large-scale statistical forecasting models reassess the unpredictability of chaotic systems. arXiv:2303.08011
  * [5] Jaeger, H. (2001). The “echo state” approach to analysing and training recurrent neural networks-with an erratum note. GMD Report, 148.
  * [6] Hara, M. & Kokubu, H. (2022). Learning Dynamics by Reservoir Computing. J. Dyn. Diff. Eq., https://link.springer.com/article/10.1007/s10884-022-10159-w
  * [7] Hirsch, M. W., Smale, S. & Devaney, R. L. (2012). Differential Equations, Dynamical Systems, and an Introduction to Chaos (third edition). Academic Press.
  * [8] Kobayashi, M. U., Nakai, K., Saiki, Y. & Tsutsumi, N. (2021). Dynamical system analysis of a data-driven model constructed by reservoir computing. Phys. Rev. E, 104, 044215.
  * [9] F. Takens, Detecting strange attractors in turbulence, in Dynamical Systems of Turburenee, Lecture Notes in Mathematics, 898, Springer-Verlag, 1981, Berlin, 366-381.
  * [10] Wold, H. O. A. (1954). A Study in the Analysis of Stationary Time Series: With an Appendix. Almqvist & Wiksell.