# A comparative study of different machine learning methods for dissipative
quantum dynamics

Luis E. Herrera Rodríguez  Department of Physics and Astronomy, University of
Delaware, Newark, DE 19716, United States  Departmento de Física, Universidad
Nacional de Colombia, Bogotá D. C., Colombia  Arif Ullah  State Key Laboratory
of Physical Chemistry of Solid Surfaces, Fujian Provincial Key Laboratory of
Theoretical and Computational Chemistry, Department of Chemistry, and College
of Chemistry and Chemical Engineering, Xiamen University, Xiamen 361005, China
Kennet J. Rueda Espinosa  Department of Physics and Astronomy, University of
Delaware, Newark, DE 19716, United States  Departmento de Física, Universidad
Nacional de Colombia, Bogotá D. C., Colombia  Pavlo O. Dral<EMAIL_ADDRESS>State Key Laboratory of Physical Chemistry of Solid Surfaces, Fujian
Provincial Key Laboratory of Theoretical and Computational Chemistry,
Department of Chemistry, and College of Chemistry and Chemical Engineering,
Xiamen University, Xiamen 361005, China  Alexei A. Kananenka<EMAIL_ADDRESS>Department of Physics and Astronomy, University of Delaware, Newark, DE 19716,
United States

###### Abstract

It has been recently shown that supervised machine learning (ML) algorithms
can accurately and efficiently predict the long-time populations dynamics of
dissipative quantum systems given only short-time population dynamics. In the
present article we benchmaked 22 ML models on their ability to predict long-
time dynamics of a two-level quantum system linearly coupled to harmonic bath.
The models include uni- and bidirectional recurrent, convolutional, and fully-
connected feed-forward artificial neural networks (ANNs) and kernel ridge
regression (KRR) with linear and most commonly used nonlinear kernels. Our
results suggest that KRR with nonlinear kernels can serve as inexpensive yet
accurate way to simulate long-time dynamics in cases where the constant length
of input trajectories is appropriate. Convolutional Gated Recurrent Unit model
is found to be the most efficient ANN model.

##  I Introduction

Simulation of the dynamics of quantum dissipative systems Weiss (2012); Breuer
and Petruccione (2002); Leggett  _et al._ (1987) is one of the most
challenging problems in physics and chemistry. Quantum dissipation arises from
the coupling of a quantum system to a thermal bath which consists of an
infinite number of degrees of freedom. This results in a time irreversible
dynamics of the system.

A multitude of numerically exact methods has been developed for quantum
dynamics simulations including the hierarchical equations of motion (HEOM),
Tanimura and Kubo (1989); Tanimura (2020) multi-configurational time-dependent
Hartree (MCTDH), Meyer  _et al._ (1990); Wang and Thoss (2003) quasi-adiabatic
propagator path integral (QUAPI), Makarov and Makri (1994) time evolving
density matrix using orthogonal polynomials algorithm (TEDOPA), Prior  _et
al._ (2010) time-dependent density matrix renormalization group (TD-DMRG), Ren
_et al._ (2018) time-dependent Davydov ansatz, Luo  _et al._ (2010)
diagrammatic quantum Monte Carlo, Cohen and Rabani (2011); Cohen  _et al._
(2015) tensor-train split-operator Fourier transform (TT-SOFT), Greene and
Batista (2017) and the stochastic equation of motion approach. Yan and Shao
(2016); Hsieh and Cao (2018a, b); Han  _et al._ (2020); Ullah  _et al._ (2020)
In practice, however, none of these methods can be used to simulate long-time
quantum dynamics of quantum systems coupled to realistic baths containing
large number of, generally, anharmonic degrees of freedom (DOF).

Methods based on Nakajima–Zwanzig generalized quantum master equation (GQME)
Nakajima (1958); Zwanzig (1960); Shi and Geva (2003); Kelly and Markland
(2013); Mulvihill and Geva (2021); Mulvihill  _et al._ (2021); Brian and Sun
(2021) including transfer tensor method (TTM), Cerrillo and Cao (2014);
Kananenka  _et al._ (2016); Buser  _et al._ (2017); Gelzinis  _et al._ (2017);
Chen  _et al._ (2020) allow to obtain long-time dynamics of the reduced
density matrix of a quantum system at a significantly lower computational cost
compared to numerically exact methods, provided the GQME memory kernels are
available. However, in general, it is difficult to numerically calculate the
exact memory kernel for GQME. TTM further reduces the computational cost
compared to the direct solution of GQME but it still requires an input set of
dynamical maps to be generated by a numerically accurate method rendering such
approaches infeasible for systems with many DOFs. Machine learning (ML) offers
an alternative route to accurate, yet greatly accelerated quantum dynamics
calculations with minimum input information required. Herrera Rodríguez and
Kananenka (2021); Ullah and Dral (2021, 2022a, 2022b)

In ML, predicting the future time-evolution of quantum mechanical observables
based on past values can be formulated as time series forecasting problem. A
time series is a set of data points recorded in consecutive time intervals.
Forecasting is a challenging part of time series data analysis. Traditional
approaches to forecasting of time series utilize moving averages. One such
approach, Autoregressive Integrated Moving Average (ARIMA) is a linear
regression-based method which, together with its variations (e.g., ARIMAX and
SARIMA), Box  _et al._ (2015) have become the standard tools for modeling
various time series problems. Ariyo  _et al._ (2014); Khashei and Bijari
(2011) These approaches perform reasonably well for short-term forecasting,
but the performance of these methods deteriorates severely for long-term
predictions. Additionally, such methods require assumptions about the
underlying data that have to be incorporated into the model.

Machine learning is artificial intelligence-based data analysis technique that
is data-driven as opposed to traditional model-driven approaches. Perhaps, the
most widely known ML tool is a feed-forward neural network (FFNN) which is an
artificial neural network (ANN) wherein connections between the nodes do not
form a cycle. Schmidhuber (2015); Goodfellow  _et al._ (2016) In this kind of
ANN the information moves only in one direction from the input nodes, through
the hidden nodes (if any) and to the output nodes without using any feedback
from the past input data. As a result, the output of any layer does not affect
the training process performed on the same layer. These types of ANNs are, in
general, not able to handle sequential inputs and all their inputs (and
outputs) must be independent of each others.

In contrast, recurrent neural networks (RNNs), Sherstinsky (2020); Hochreiter
and Schmidhuber (1997); Goodfellow  _et al._ (2016) in general, are designed
with built-in gates that function to store the previous inputs and leverage
sequential information. RNNs utilizes a loop in the network to preserve some
information and thus functions like a memory. This gives such feedback-based
models the ability to learn from past data. RNNs are known to outperform
ARIMA-based models in the long-time forecasting problems. Ye  _et al._ (2012);
Siami-Namini  _et al._ (2018, 2019)

RNNs have been successfully applied to time series classification Malhotra
_et al._ (2017); Fawaz  _et al._ (2019); Kashiparekh  _et al._ (2019) and
forecasting across many domains including financial data predictions, Moghar
and Hamiche (2020); Fischer and Krauss (2018); Kim and Moon (2019); Kumar  _et
al._ (2021) speech Graves (2012a); Graves and Schmidhuber (2005); Graves
(2012b); Graves  _et al._ (2013); Hannun  _et al._ (2014); Xiong  _et al._
(2016) and handwriting recognition, Graves  _et al._ (2007); Graves and
Schmidhuber (2009) natural language processing, Tang  _et al._ (2016); Yin
_et al._ (2017) machine translation, Sutskever  _et al._ (2014); Bahdanau  _et
al._ (2014); Shido  _et al._ (2019) healthcare, Choi  _et al._ (2017); Lynn
_et al._ (2019); Gupta  _et al._ (2021); Ma  _et al._ (2017) traffic speed
prediction, van Lint  _et al._ (2005); Zhao  _et al._ (2017); Cui  _et al._
(2019) music, Eck and Schmidhuber (2002); Boulanger-Lewandowski  _et al._
(2012) video, Srivastava  _et al._ (2015) meteorology, Habi and Messer (2019)
molecular drug discovery, Gupta  _et al._ (2018); Segler  _et al._ (2018)
demand forecasting, Wei and Yan (2019); Abbasimehr  _et al._ (2020) gaming,
Zhang  _et al._ (2020) remote sensing Lyu  _et al._ (2016); Ienco  _et al._
(2017) computer code generation Bhoopchand  _et al._ (2016) and others.

There are several variations of RNN-based models differing in their
capabilities to remember input data. The vanilla RNN, Goodfellow  _et al._
(2016); Pineda (1987) long short-term memory (LSTM), Hochreiter and
Schmidhuber (1997); Graves and Schmidhuber (2005); Gers  _et al._ (2000,
1999); Goodfellow  _et al._ (2016) and gated recurrent unit (GRU) Cho  _et
al._ (2014); Goodfellow  _et al._ (2016) are the most commonly used types of
RNNs. The vanilla RNN is the first model of recurrent ANNs that was
introduced. Pineda (1987) Its ability to learn long-term dependencies is
limited due to vanishing and exploding gradient problems. LSTM-based models
are extensions of RNNs with significantly improved performance on long-term
predictions. LSTM uses a set of gates to learn what information worth to
remember. GRU is another gated variant of RNN model developed by Chung et al.
Chung  _et al._ (2016) who showed that GRU outperforms the LSTM in some tasks.
Josefowicz et al. Jozefowicz  _et al._ (2015) also reported the superior
performance of the GRU-based models, but noticed that the performance of LSTM
can nearly match that of the GRU if proper initialization is performed.

Bidirectional RNNs (BRNN) Schuster and Paliwal (1997) are another extension of
the standard unidirectional RNNs in which two RNNs are applied to the input
data. Firstly, an RNN is applied on the input sequence which is then followed
by the application of RNN on the reversed input sequence. This typically
improves the accuracy of the model Baldi  _et al._ (1999); Wei and Yan (2019)
at a cost of slower training times. Siami-Namini  _et al._ (2019); Wei and Yan
(2019) Several variants of bidirectional RNNs exist differing in the type of
the underlying RNN cell. Bidirectional LSTMs (BLSTM) Graves  _et al._ (2013)
and bidirectional GRUs (BGRU) Lynn  _et al._ (2019) are being explored in
various tasks. Graves  _et al._ (2013); Tang  _et al._ (2016); Ma  _et al._
(2017); Kim and Moon (2019); Wei and Yan (2019); Zhang  _et al._ (2020) In
some applications, Siami-Namini  _et al._ (2019); Cui  _et al._ (2019) such as
speech recognition, the better performance of BLSTM compared to the regular
unidirectional LSTM has been reported Graves  _et al._ (2013) and is not
surprising given the nature of the task (text parsing and prediction of next
words in the input sentence). In general, however, it is not clear what
problems benefit from the bi-directional training.

Convolutional neural networks (CNNs) LeCun  _et al._ (1989) are a special kind
of ANNs designed for processing data that has a known grid-like topology.
Goodfellow  _et al._ (2016) For example, time-series data can be thought of as
a one-dimensional (1D) grid taking samples at regular time intervals. One-
dimensional CNNs (1D CNNs) have achieved promising results in time-series
classification tasks Wang  _et al._ (2016); Cui  _et al._ (2016); Serrà  _et
al._ (2018); Zheng  _et al._ (2014, 2016); Interdonato  _et al._ (2018);
Kashiparekh  _et al._ (2019); Ismail Fawaz  _et al._ (2020); Tang  _et al._
(2021) in many domains including healthcare, Rajpurkar  _et al._ (2017); Roy
_et al._ (2018); Schirrmeister  _et al._ (2018) speech recognition, Sercu  _et
al._ (2016); Xiong  _et al._ (2016) music classification, Choi  _et al._
(2016); You  _et al._ (2018) natural language processing Yin  _et al._ (2017)
and others.

CNNs have also been combined with RNNs. Such hybrid models are called
convolutional recurrent networks (CRNNs). A CRNN is a deep ANN that contains a
CNN layer(s) followed by an RNN layer(s). Such architectures possess several
advantages. Firstly, 1D CNN layers learn to sub-sample the data and reduce the
input vector that is passed to an RNN layer(s). This is important because GRU
and LSTM layers are computationally expensive and replacing some of them with
convolutional layer(s) improves the computational scaling of the algorithm.
Second, 1D CNN layers extract local information from neighboring time points
and pass already detected temporal dependencies further down to RNN layers.
CRNNs are being actively explored in various tasks. Karim  _et al._ (2018);
Roy  _et al._ (2018); You  _et al._ (2018); Choi  _et al._ (2016); You  _et
al._ (2018); Sheykhivand  _et al._ (2020) It has been reported that CRNNs can
outperform CNNs in some tasks. Choi  _et al._ (2016) CNNs were combined with
BLSTMs as well. Eapen  _et al._ (2019) Note that these methods are different
from convolutional LSTM models Shi  _et al._ (2015) in which input
transformations and recurrent transformations are both convolutional.

Kernel methods represent another class of ML methods that are applied to time-
series analysis. Müller  _et al._ (1997); Sapankevych and Sankar (2009);
Haworth  _et al._ (2014) Such methods employ a function (kernel) that maps
input data into a high dimensional space and then perform a linear regression
in that space. Kernel Ridge Regression (KRR) and Support Vector Regression
(SVR) are examples of such algorithms. A crucial aspect of applying kernel
methods to time series data is to find appropriate kernels to distinguish
between time series. A simple way is to treat time series as static vectors,
essentially as they are treated in a feed-forward ANN, ignoring the time
dependence. In such cases standard kernels such as Matérn Gneiting  _et al._
(2010) or Gaussian radial basis function (RBF) Gneiting  _et al._ (2010) can
be used. However, such methods are limited to input time-sequences of equal
length, again similar to ANNs, yet many applications involve time sequences of
varying length. To overcome this limitation, kernel functions such as
autoregressive Cuturi and Doucet (2011) kernel have been developed.

Recently, many ML models have been applied to simulate the dynamics of quantum
systems. Herrera Rodríguez and Kananenka (2021); Ullah and Dral (2021, 2022b,
2022a); Akimov (2021); Secor  _et al._ (2021); Yang  _et al._ (2020);
Bandyopadhyay  _et al._ (2018); Banchi  _et al._ (2018); Wu  _et al._ (2021);
Lin  _et al._ (2021); Choi  _et al._ (2022); Tsai  _et al._ (2022); Tang  _et
al._ (2022) We note that ML can also be applied to quantum dynamics in a
different context—namely as surrogate models for quantum chemical properties
such as potential energies and forces in different electronic states as well
as couplings between states eliminating the need for expensive (excited-state)
electronic structure calculations. Dral and Barbatti (2021); Westermayr and
Marquetand (2020, 2021) Here we apply ML to propagate a quantum system
assuming potential energies are readily available. ML models are attractive
because of their very low computational cost and favorable scaling with
respect to the size of a quatum system and bath. RNNs were used to simulate
dynamics of the spin-boson model Hamiltonian and Landau–Zener transitions,
Bandyopadhyay  _et al._ (2018); Yang  _et al._ (2020) and to learn the
convolutionless master equation. Banchi  _et al._ (2018) Rodriguez et al.
Herrera Rodríguez and Kananenka (2021) used CNNs to accurately model long-time
dynamics of the spin-boson system based on short-time dynamical data. Later
Ullah et al. Ullah and Dral (2021) illustrated that KRR methods can also
predict long-time dynamics of the spin-boson model very accurately. Recently,
CNNs is used to study the excitation energy transfer in Fenna–Matthews–Olson
light-harvesting complex. Ullah and Dral (2022a, b) Wu et al. Wu  _et al._
(2021) used a hybrid CNN/LSTM network to predict long-time semiclassical and
mixed quantum-classical dynamics of the spin-boson model. Lin et al. Lin  _et
al._ (2021, 2022) trained a multi-layer LSTM model to simulate the long-time
dynamics of spin-boson model and used bootstrap method to estimate the
confidence interval. We note that ML methods based on FFNNs have also been
recently applied to model quantum dynamics. Akimov (2021); Secor  _et al._
(2021) The recent upsurge of applications of ML methods to dissipative quantum
dynamics calls for a systematic benchmark of such methods.

In this article, we present a comprehensive comparison of 22 ML models for
predicting the long-time dynamics of an open quantum system given the short-
time evolution data. We consider all three most used types of unidirectional
RNNs including the vanilla RNN, GRU, and LSTM, the corresponding bidirectional
RNNs (BRNN, BGRU, BLSTM), 1D CNNs, CRNNs, as well as KRR with 8 different
kernel functions. We test the ability of these ML methods to predict the
donor-acceptor population difference of a spin-boson model in symmetric and
asymmetric regimes over a broad range of temperatures, reorganization
energies, and bath relaxation timescales.

##  II Theory

###  II.1 Model system

We choose to test ML methods on the spin-boson model which has become a
paradigmatic model system in the study of open quantum systems due to the
richness of its physics. Leggett  _et al._ (1987) The spin-boson model was
exploited in a wide-range of applications in quantum computing, Makhlin  _et
al._ (2001) quantum phase transitions, Alvermann and Fehske (2009); Winter
_et al._ (2009) electron transfer in biological systems, Garg  _et al._ (1985)
and others. The spin-boson model describes a two-level quantum subsystem
linearly coupled to a heat-bath environment. The bath is modeled as an
ensemble of independent harmonic oscillators. The total Hamiltonian in the
subsystem’s basis $\\{|+\rangle,|-\rangle\\}$ is given by ($\hbar=1$)

$\hat{H}=\frac{\epsilon}{2}\hat{\sigma}_{z}+\frac{\Delta}{2}\hat{\sigma}_{x}+\hat{\sigma}_{z}\sum_{\alpha}g_{\alpha}\left(b_{\alpha}^{\dagger}+b_{\alpha}\right)+\sum_{\alpha}\omega_{\alpha}b_{\alpha}^{\dagger}b_{\alpha},$
(1)

where $b^{\dagger}_{\alpha}$ ($b_{\alpha}$) is the bosonic creation
(annihilation) operator of the $\alpha$th mode with frequency
$\omega_{\alpha}$, $\hat{\sigma}_{z}=|+\rangle\langle+|-|-\rangle\langle-|$
and $\hat{\sigma}_{x}=|+\rangle\langle-|+|-\rangle\langle+|$ are the Pauli
operators, $\epsilon$ is the energetic bias, $\Delta$ is the tunneling matrix
element, and $g_{\alpha}$ are the coupling coefficients.

The impact of the bath is completely determined by the spectral density
$J(\omega)=\pi\sum_{\alpha}g_{\alpha}^{2}\delta(\omega_{\alpha}-\omega)$
which, in this work, is choosen to be of the Debye form (Ohmic spectral
density with the Drude–Lorentz cut-off) Wang  _et al._ (1999)

$J(\omega)=2\lambda\frac{\omega\omega_{c}}{\omega^{2}+\omega_{c}^{2}},$ (2)

where $\lambda$ is the bath reorganization energy which controls the strength
of the coupling between system and the bath, and the cutoff frequency
$\omega_{c}$ which sets the primary timescale for the bath evolution
$\tau_{c}=(\omega_{c})^{-1}$.

We consider the time evolution of the expectation value of $\hat{\sigma}_{z}$
Pauli operator

$\langle{\sigma}_{z}(t)\rangle=\mathrm{Tr}_{\mathrm{s}}\left[\hat{\sigma}_{z}\hat{\rho}_{\mathrm{s}}(t)\right],$
(3)

which is often referred to as the population difference
$\langle{\sigma}_{z}(t)\rangle=p_{+}(t)-p_{-}(0)$, where
$p_{\pm}(t)=\mathrm{Tr}_{\mathrm{s}}\left[|\pm\rangle\langle\pm|\hat{\rho}_{s}(t)\right]$.
In Eq. (3) the trace is taken over the system degrees of freedom as denoted by
“s”, and $\hat{\rho}_{\mathrm{s}}$ is the system’s reduced density operator

$\hat{\rho}_{\mathrm{s}}(t)=\mathrm{Tr}_{\mathrm{b}}\left[e^{-i\hat{H}t}\hat{\rho}(0)e^{i\hat{H}t}\right],$
(4)

where $\hat{\rho}(0)$ is the total system plus bath density operator and the
trace is taken over the bath degrees of freedom. The initial state of the
total system is assumed to be a product state of the following form

$\hat{\rho}(0)=\hat{\rho}_{\mathrm{s}}(0)\otimes\frac{e^{-\beta\hat{H}_{\mathrm{b}}}}{Z_{\mathrm{b}}},$
(5)

where
$Z_{\mathrm{b}}=\mathrm{Tr}_{\mathrm{b}}\left[e^{-\beta\hat{H}_{\mathrm{b}}}\right]$
is the bath partition function, $\beta=(k_{\mathrm{B}}T)^{-1}$ is the inverse
temperature, and $k_{\mathrm{B}}$ is the Boltzmann constant. The initial
density operator of the system is chosen to be
$\hat{\rho}_{\mathrm{s}}(0)=|+\rangle\langle+|$. These conditions correspond
to situations where the initial preparation of the subsystem occurs quickly on
the timescale of the bath relaxation.

###  II.2 Machine learning models

In this section we provide a detailed description of all ML models used in the
present study. We specialize our discussion to modeling time-series data with
the input sequence denoted as $\mathbf{x}=\left(x^{(1)},\ldots,x^{(T)}\right)$
where $T$ is the length of the time series. Each element $x^{(t)}$ of
$\mathbf{x}$ can be a real-valued vector itself, $x^{(t)}\in\mathbb{R}^{n}$.
In this work, the dimension $n$ of each element of an input sequence
$\mathbf{x}^{(t)}$ is 1. Consider a data set
$\mathcal{D}=\\{(\mathbf{x}_{i},\mathbf{y}_{i})\\}_{i=1}^{N}$ containing $N$
time series $\mathbf{x}_{i}$ and their associated labels $\mathbf{y}_{i}$. In
time-series forecasting problems, labels can describe the future states of the
input sequence $\mathbf{x}_{i}$ as denoted by
$\mathbf{y}_{i}=\left(x_{i}^{(T+1)},\ldots,x_{i}^{(T+m)}\right)$ which
corresponds to $m$ next elements of the sequence $\mathbf{x}_{i}$. Time-series
forecasting in a supervised learning framework amounts to (i) training ML
models on the subset of $\mathcal{D}$ called a training set and (ii) using
trained ML models to make a prediction $\hat{\mathbf{y}}_{i}$ for a given test
time series $\mathbf{x}_{i}$. In this work we test the ability of ML models to
predict a single real-valued scalar quantity, the population difference of the
spin-boson model, $\langle{\hat{\sigma}}_{z}(t)\rangle$ for a single time step
i.e., $m=1$. Extension to multivariate time series data ($n>1$) and multi-step
outputs ($m>1$) is possible.Ullah and Dral (2022b)

####  II.2.1  Feedforward neural networks

Feedforward neural networks (FFNNs) or multiplayer perceptrons are the most
used type of artificial neural networks. FFNN approximate a function
$g(\mathbf{x})$ by defining a mapping
$\mathbf{y}=\mathcal{F}(\mathbf{x};\bm{\theta})$ and determining the value of
parameters $\bm{\theta}$ that best approximate $g(\mathbf{x})$. These models
are called feedforward because information flows through the function being
evaluated from an input $\mathbf{x}$ through intermediate steps to the output
$\hat{\mathbf{y}}$. There are no feedback connections between the output and
the input of the model. Goodfellow  _et al._ (2016)

A crucial element of the success of ANNs is the use of deep architectures.
Deep ANNs are created by stacking multiple layers on top of each other, with
the output of one layer forming the input sequence for the next layer. The
first layer is called an input layer, the last layer is called the output
layer; all layers in between are called hidden layers.

Deep FFNNs are compositions of several functions:
$\hat{\mathbf{y}}=\mathcal{F}^{(L)}\large(\ldots\mathcal{F}^{(2)}(\mathcal{F}^{(1)}\left(\mathbf{x};\bm{\theta}^{(1)})\right);$
$\bm{\theta}^{(2)});\bm{\theta}^{(L)}\large)$ where each function
$\mathcal{F}^{(k)}$ depends on its own set of parameters $\bm{\theta}^{(k)}$.
Function $\mathcal{F}^{(l)}$ connects $l$th and $(l+1)$th layers of the
network. It takes an input $\mathbf{x}^{(l)}\in\mathbb{R}^{k_{l}}$ and
generates the output according to

$\mathbf{x}^{(l+1)}=\mathcal{F}^{(l)}\left(\mathbf{x}^{(l)};\bm{\theta}^{(l)}\right)=f^{(l)}\left(\mathbf{a}^{(l)}\right),$
(6)

where $f^{(l)}:\mathbb{R}\to\mathbb{R}$ is the $l$th layer activation function
which is applied elementwise. As seen from Eq. (6) each layer of the network
is vector valued. Each vector element $a_{j}^{(l)}$ is called a neuron. Its
value is calculated from layer’s input and parameters

$a_{j}^{(l)}=\sum_{i=1}^{k_{l}}w_{ij}^{(l)}x_{i}^{(l)}+b_{j}^{(l)},$ (7)

where $\mathbf{w}^{(l)}\in\mathbb{R}^{k_{l}\times k_{l+1}}$ and
$\mathbf{b}^{(l)}\in\mathbb{R}^{k_{l}}$ are called the weights and the biases
of the $l$th layer, respectively, that together constitute the set of
trainable parameters of the $l$th layer
$\bm{\theta}^{(l)}=\\{\mathbf{w}^{(l)},\mathbf{b}^{(l)}\\}$. The output
$\mathbf{x}^{(l+1)}$ forms an input into $(l+1)$th layer. The total number of
trainable parameters, biases and weights, of a single FFNN layer is
$k_{l}(k_{l+1}+1)$. It should be emphasized that FFNN models require input
time sequences of a fixed length, $T$. This, in particular, requires an a
priori knowledge of the memory effects in a quantum system under study which
is non-trivial task. Artificial neural networks based on recurrent layers,
discussed below, do not impose such a restriction.

The strategy of deep learning is to find the set of model parameters
$\\{\bm{\theta}^{(1)},\ldots,\bm{\theta}^{(L)}\\}$ that best approximates the
target function $g(\mathbf{x})$. The model parameters are adjusted during the
training process which is based on backpropagation algorithm. Goodfellow  _et
al._ (2016)

Typically the activation function of the hidden layers is chosen to be
nonlinear, e.g., logistic sigmoid function $\sigma(z)=(1+e^{-z})^{-1}$.
According to the universal approximation theorem Hornik  _et al._ (1989);
Cybenko (1989); Leshno  _et al._ (1993) an FFNN with a linear output
activation function and at least one hidden layer with a nonlinear activation
function can approximate any Borel measurable function (continuous on the
closed and bonded subset of the real coordinate space) from one finite-
dimensional space to another with any desired nonzero amount of error,
provided the network contains enough neurons. The theorem gurantees that
regardless of the target function, a single hidden layer FFNNs with many
neurons will be able to represent this function with any degree of accuracy.
However, in general, there is no guarantee that the training algorithm can do
so. Goodfellow  _et al._ (2016)

####  II.2.2 1D convolutional neural networks

Convolutional neural networks (CNNs) is a type of ANNs that can be applied to
time-series data. LeCun  _et al._ (1989); Goodfellow  _et al._ (2016) CNNs are
based on a mathematical operation called convolution. Let
$\left(\mathbf{g}\ast\mathbf{w}\right)$ be the result of a 1D discrete
convolution and the $i$th element of the result is given by

$\left(\mathbf{g}\ast\mathbf{w}\right)_{i}=\sum_{j=1}^{k}g_{i+j-1}w_{k+1-j},$
(8)

where $\mathbf{w}\in\mathbb{R}^{k}$ is referred to as kernel, or weight
vector, and $\mathbf{g}$ is a time series. CNNs exploit two key ideas: sparse
connectivity and parameter sharing. Goodfellow  _et al._ (2016) The former is
accomplished by making the kernel size smaller than the size of the input
which allows to detect, for the time-series data, short-time correlations and
reduces the number of parameters to be stored in memory compared to FFNNs. The
latter amounts to using the same kernel parameters for all positions in the
input which further reduces the storage requirements. In the case of
convolution, the parameter sharing leads to equivariance—a property that
causes the output to change in the same way the input changes. Specifically,
in the case of time-series data, the convolution captures a timeline for
different features to appear in the input. Goodfellow  _et al._ (2016)

In practical implementations of CNNs, the convolution operation given in Eq.
(8) is often replaced by the so-called cross-correlation operation Goodfellow
_et al._ (2016) which we denote by a “$\star$”

$\left(\mathbf{g}\star\mathbf{w}\right)_{i}=\sum_{j=1}^{k}g_{i+j-1}w_{j}.$ (9)

It should be noted that Eq. (9) is given for the step size of the kernel, as
it is applied across the sequence, equal to 1. This step size is called stride
$s(s\in\mathbb{Z}^{+})$. In general, for the stride $s\geq 1$, Eq. (9) is
modified to

$\left(\mathbf{g}\star\mathbf{w}\right)_{i}=\sum_{j=1}^{k}g_{(i-1)s+j}w_{j}.$
(10)

Let’s consider a two-layer 1D CNN architecture. An input sequence $\mathbf{x}$
is processed through $K_{1}$ 1D kernels of the first layer whose weights are
denoted as $\mathbf{w}^{(1)}\in\mathbb{R}^{k_{1}}$. The result of this
operation is a tensor $\mathbf{Z}^{(1)}$ whose elements are given by

$Z_{ij}^{(1)}=f^{(1)}\left[\left(\mathbf{x}\star\mathbf{w}^{(1)}_{j}\right)_{i}+b_{j}^{(1)}\right],$
(11)

where $f^{(1)}$ is the activation function,
$\mathbf{b}^{(1)}\in\mathbb{R}^{K_{1}}$ are the bias parameters and
$\mathbf{w}_{j}^{(1)},j=1,\ldots,K_{1}$ is the $j$th kernel of the first
layer. The output tensor $\mathbf{Z}^{(1)}$ has dimensions of $M_{1}\times
K_{1}$ with $M_{1}=\lfloor(T-k_{1}+2p_{1})/s_{1}\rfloor+1$, where $s_{1}$ is
the stride, $p_{1}$ is the amount of zero padding, and $\lfloor
x\rfloor=\max\\{m\in\mathbb{Z}|m\leq x\\}$ is the floor function. Zero padding
parameter $p_{1}$ denotes the number of zeros to be added to the start and the
end of the input sequence. For example, for $p=0$ the input sequence is not
padded with zeros and, consequently, the kernel is allowed to visit only
positions where it is contained entirely within the input sequence. This type
of zero padding is called a valid padding.

Because convolution of input data with a single kernel can extract only one
kind of features, albeit at many locations, in practice, many filters are
used. Each element of a kernel is a trainable parameter and the same
parameters of each kernel are shared across the entire sequence. Therefore,
the total number of trainable parameters of the first 1D CNN layer is
$K_{1}(k_{1}+1)$.

The output of the first layer $\mathbf{Z}^{(1)}$ is then processed through
$K_{1}\times K_{2}$ 1D kernels whose weights are denoted as
$\mathbf{w}^{(2)}\in\mathbb{R}^{k_{2}}$. The output of the second layer is
tensor $\mathbf{Z}^{(2)}$ whose elements are given by

$Z_{ij}^{(2)}=f^{(2)}\left[\sum_{u=1}^{K_{1}}\left(\mathbf{Z}^{(1)}_{u}\star\mathbf{w}_{uj}^{(2)}\right)_{i}+b_{j}^{(2)}\right],$
(12)

where $f^{(2)}$ is the activation function of the second layer and
$\mathbf{b}^{(2)}$ are the bias parameters of the second layer. In Eq. (12)
$\mathbf{w}_{uj}^{(2)}$ is to be understood as the the $u$th kernel (weight
vector) applied to the output of $j$th kernel (also called channel) of the
first layer. The dimensions of $\mathbf{Z}^{(2)}$ are $M_{2}\times K_{2}$
where $M_{2}=\lfloor(M_{1}-k_{2}+2p_{2})/s_{2}\rfloor+1$ with $p_{2}$ and
$s_{2}$ being the zero-padding and the stride of the second 1D CNN layer,
respectively. It is worth emphasizing that a separate set of kernels is
applied to the output of each kernel of the preceding (in this case first)
layer. Thus, the total number of trainable parameters of the second 1D CNN
layer is $K_{2}(K_{1}k_{2}+1)$.

A pooling layer is commonly found in CNN architectures. Its function is to
subsample (shrink) the input. Pooling function replaces the output of a CNN
layer at a certain location with a neighborhood-dependent information. For
example, the maximum pooling (MaxPooling) operation Zhou and Chellappa (1988)
outputs the maximum value in a neighborhood of a specified size. This helps to
make representation approximately invariant to small translations of the input
and improves the computational efficiency in a deep CNN by reducing the number
of trainable parameters of the next convolutional layer. A 1D MaxPooling layer
with the kernel size $k_{p}$, stride $s_{p}$, operates independently on every
depth slice of the input $\mathbf{Z}^{\mathrm{in}}$ and resizes it by taking
only the maximum value

$Z_{ij}^{\mathrm{out}}=\max\left(Z_{(i-1)s_{p}+1,j}^{\mathrm{in}},\ldots,Z_{(i-1)s_{p}+k_{p},j}^{\mathrm{in}}\right),$
(13)

where $i=1,\ldots,\lfloor(M_{2}-k_{p})/s_{p}\rfloor+1$. Note that Eq. (13) is
written, for simplicity, for the zero padding but other kinds of padding can
be used as well. MaxPool layers do not employ any trainable parameters.

####  II.2.3 Recurrent neural networks

Recurrent neural networks (RNNs) are a family of neural networks specifically
designed for processing sequential data. Rumelhart  _et al._ (1986) Similarly
to FFNNs, RNNs have the universal approximation ability. Schäfer and
Zimmermann (2006); Hornik  _et al._ (1989) Similarly to CNNs, RNNs are based
on the idea of parameter sharing but, unlike CNNs, RNNs share parameters
through recursion. Formally, given a sequence $\mathbf{x}$ at each time step,
an RNN updates its hidden state
$\mathbf{H}=\left[\mathbf{h}^{(1)},\ldots,\mathbf{h}^{(T)}\right]$ recursively
based on the current input $x^{(j)}$ and the previous hidden state
$\mathbf{h}^{(j-1)}$ as follows

$\mathbf{h}^{(j)}=\mathcal{R}(\mathbf{h}^{(j-1)},x^{(j)}),$ (14)

where $\mathcal{R}$ is a nonlinear function. Each element of the hidden state
vector, $\mathbf{h}^{(j)}$ is, in general, a vector itself. Recursive
application of Eq. (14) results in the sharing of parameters across an ANN
architecture. Training of such ANN amounts to selectively emphasizing some
aspects of the past sequence inputs that deemed more important than others.

All RNN models used in this work contain two stacked recurrent layers and are
schematically illustrated in Fig. 1a. The number of RNN cells in the first
layer equals to the length of an input vector such that each cell processes a
single time step $t=1,\ldots,T$ and updates the corresponding hidden state
vector $\mathbf{h}^{(t)}\in\mathbb{R}^{k_{1}}$. The size of the hidden state
vector $k_{1}$ is a hyperparameter that needs to be optimized. Note that
various software packages use different terminology. In this work RNN models
are built with the Keras Chollet  _et al._ (2015) software whose argument
units corresponds to the size of the hidden state vector. An RNN layer outputs
a tensor, $\mathbf{H}\in\mathbb{R}^{T\times k_{1}}$, containing all hidden
states and it serves as an input to the following RNN layer.

Each RNN cell of a layer processes the corresponding slice of the preceding
layer’s output $\mathbf{H}_{Lj}$, where $j\in[1,k_{1}]$ and $L$ is the number
of cells of the second layer, and transforms it into the state vector
$\mathbf{h}_{2}^{(t)}\in\mathbb{R}^{k_{2}}$, where $k_{2}$ is the number of
units of the second RNN layer. The last RNN layer can either return the hidden
state vector at only the final time step $T$ or the entire sequence of hidden
states.

Various RNN architectures differ by the function $\mathcal{R}$ in Eq. (14). A
simple RNN cell, also known as vanilla RNN is illustrated in Fig. 1b with the
hyperbolic tangent activation function. It receives an input $x^{(t)}$
corresponding to the time step $t$ and uses the state vector from the
chronologically previous RNN cell $\mathbf{h}^{(t-1)}\in\mathbb{R}^{k}$ to
generate the hidden state $\mathbf{h}^{(t)}$ according to the following update
equation

$\mathbf{h}^{(t)}=\tanh\left(\mathbf{b}+\mathbf{w}_{h}\cdot\mathbf{h}^{(t-1)}+\mathbf{w}_{x}\cdot\mathbf{x}^{(t)}\right),$
(15)

where $\mathbf{b}_{h}\in\mathbb{R}^{k}$ are the bias parameters,
$\mathbf{w}_{h}\in\mathbb{R}^{k\times k}$,
$\mathbf{w}_{x}\in\mathbb{R}^{k\times m}$ are the coefficient matrices, and
the tanh activation function is applied element-wise. The function
$\mathcal{R}(\mathbf{h}^{(j-1)},x^{(j)})$ in case of the vanilla RNN amounts
to iterative application of Eq. (15) to an input sequence $\mathbf{x}$ from
$t=1$ to $t=T$.

Figure 1: (a) Basic structure of a two-layer unidirectional recurrent neural
network. Blue rounded rectangles represent cells and $x_{t}$ is the $t$th
element of an input sequence $\mathbf{x}$. Hidden state vector
$\mathbf{h}^{(t)}$ is passed between the cells within a layer as well as
between the layers; (b) the vanilla RNN cell with tanh activation function;
(c) LSTM cell; (d) GRU cell. Green and yellow circles denote hyperbolic
tangent and sigmoid activation functions correspondingly. Dark blue circles
denote element-wise arithmetic operations shown on the circle. See text for
details.

The hidden state $\mathbf{h}^{(t)}$ is then passed to the next RNN cell.
Weights $\mathbf{w}_{h},\mathbf{w}_{x}$ and biases $\mathbf{b}$ in Eq. (15)
are updated iteratively via the back-propagation algorithm. In a multilayer
RNN architecture all hidden state vectors $\mathbf{H}$ are passed to the next
layer. For RNNs whose output is used directly to predict the next value of a
sequence the hidden state corresponding to $t=T$ is the desired predicted
value of the next time step $\hat{\mathbf{y}}=\mathbf{h}^{(T)}$.

The total number of trainable parameters of the vanilla RNN cell is
$k(k+n+1)$, where $k$ is the size of the hidden state vector, which is the
same for all RNN cells within a layer, and $n$ is the dimension of each
element of an input sequence $\mathbf{x}^{(t)}$ which, in this work is 1 for
the first RNN layer and $k_{1}$ for the second RNN layer. Note that if RNN is
not the first layer, as e.g., in convolutional recurrent neural networks
discussed below, $n$ would be different from 1. It should also be noted that
the total number of trainable parameters of the whole RNN layer comprised of
any number of the vanilla RNN cells is the same because
$\mathbf{b},\mathbf{w}_{h}$, and $\mathbf{w}_{x}$ are shared across all RNN
cells within the given layer. The total number of trainable parameters in a
two-layer vanilla RNN architecture is $k_{2}(k_{1}+k_{2}+1)+k_{1}(k_{1}+1+1)$.

In spite of being simple and powerful ANN model, the vanilla RNN described
above is plagued by the vanishing gradient and exploding gradient problems
when trained via backpropagation. Goodfellow  _et al._ (2016); Hochreiter
(1991); Bengio  _et al._ (1994, 1993) The exploding gradients problem refers
to the large increase in the norm of the gradient during training. The
vanishing gradients problem refers to the opposite behavior, when long term
components decay exponentially fast to zero norm, limiting the model’s ability
to learn long-range dependencies. The exploding gradient problem can be easily
addressed by gradient clipping. Pascanu  _et al._ (2013) The vanishing
gradient problem is, however, much more difficult to address.

####  II.2.4 Long short-term memory networks

In general, learning long-term dependencies is one of the most important
problems in deep learning. Goodfellow  _et al._ (2016) The difficulty in
dealing with long-memory sequences in RNNs arises from the exponentially
smaller weights given to long-term correlations compared to short-term ones.
Goodfellow  _et al._ (2016); Hochreiter (1991); Bengio  _et al._ (1994, 1993)
This problem is particular to simple RNN cell described above. It was shown
that very deep FFNNs can avoid the vanishing and exploding gradient problems
Sussillo and Abbott (2014) but, in order to store memories, RNNs must enter a
region of parameter space where gradients vanish. Bengio  _et al._ (1993,
1994)

To overcome the vanishing gradient problem various gated RNN architectures
such as long short-term memory (LSTM) were developed. Gated RNNs are based on
the idea of creating paths through time such that the derivatives are neither
vanish nor explode. The Long Short-Term Memory (LSTM) Hochreiter and
Schmidhuber (1997); Graves and Schmidhuber (2005); Gers  _et al._ (2000, 1999)
model was developed by Hochreiter et al. Hochreiter and Schmidhuber (1997) and
is based on the idea of using self-loops to produce paths where the gradient
can flow for long duration. This is achieved by adding the cell state denoted
by $\mathbf{C}^{(t)}\in\mathbb{R}^{k}$ and data-dependent gates, that control
the flow of information, to the standard RNN architecture. Hochreiter and
Schmidhuber (1997); Gers  _et al._ (2000) Both hidden state and cell state
control the memory of the network. The cell state carries relevant information
throughout the processing of the sequence such that even information from the
earlier time steps can make its way to later time steps, reducing the effects
of short-term memory. The information is added or removed to the cell state
via gates. The gates are different simplest FFNNs that decide which
information is allowed on the cell state. Thus, the gates can learn what
information is relevant to keep or forget during training.

Although a number of variants of the LSTM cell have been produced, a large-
scale analysis shows that none of them outperforms the standard LSTM
architecture. Hochreiter and Schmidhuber (1997); Graves and Schmidhuber
(2005); Gers  _et al._ (2000) The block diagram of the LSTM cell is shown in
Fig. 1c. It contains three gates: the forget gate, the input gate, and the
output gate. All gates have a sigmoid [logistic $\sigma(z)=(1+e^{-z})^{-1}$],
nonlinear activation function (yellow circles). The t$th$ element(s) of the
input sequence $\mathbf{x}^{(t)}$ enters the cell and is concatenated with the
hidden state vector from the chronologically previous cell
$\mathbf{h}^{(t-1)}$. The total vector is passed through the forget gate

$\mathbf{F}^{(t)}=\sigma\left(\mathbf{b}^{F}+\mathbf{W}_{h}^{F}\cdot\mathbf{h}^{(t-1)}+\mathbf{W}_{x}^{F}\cdot\mathbf{x}^{(t)}\right),$
(16)

where $\mathbf{b}^{F}\in\mathbb{R}^{k}$ are the biases and
$\mathbf{W}_{h}^{F}\in\mathbb{R}^{k\times k}$,
$\mathbf{W}_{x}^{F}\in\mathbb{R}^{k\times m}$, are the corresponding weights.
The forget gate was the crucial addition to the original LSTM cell extending
the length of sequences that can be processed. Gers  _et al._ (2000) The
forget gate learns to reset memory blocks once their contents are out of date
and hence useless. The output vector $\mathbf{F}^{(t)}\in\mathbb{R}^{k}$
contains the values between 1 and 0 emphasizing or diminishing the importance
of the elements of the hidden state vector, correspondingly. In two separate
branches, the current input vector and the previous hidden state vector are
processed through the input gates: the external input gate

$\mathbf{I}^{(t)}=\sigma\left(\mathbf{b}^{I}+\mathbf{W}_{h}^{I}\cdot\mathbf{h}^{(t-1)}+\mathbf{W}_{x}^{I}\cdot\mathbf{x}^{(t)}\right),$
(17)

and the new candidate gate

$\mathbf{G}^{(t)}=\tanh\left(\mathbf{b}^{G}+\mathbf{W}_{h}^{G}\cdot\mathbf{h}^{(t-1)}+\mathbf{W}_{x}^{G}\cdot\mathbf{x}^{(t)}\right),$
(18)

where $\mathbf{b}^{I}\in\mathbb{R}^{k}$, $\mathbf{b}^{G}\in\mathbb{R}^{k}$,
$\mathbf{W}_{h}^{I}\in\mathbb{R}^{k\times k}$,
$\mathbf{W}_{h}^{G}\in\mathbb{R}^{k\times k}$,
$\mathbf{W}_{x}^{I}\in\mathbb{R}^{k\times m}$, and
$\mathbf{W}_{x}^{G}\in\mathbb{R}^{k\times m}$ are the corresponding biases and
weights of the input and new candidate gates. The external output gate uses
the sigmoid activation function to obtain a gating value between 0 and 1 which
is then element-wise multiplied by the corresponding element of the new
candidate vector $\mathbf{G}^{(t)}$ effectively deciding which elements of the
input vector are the most important. Given $\mathbf{F}^{(t)}$,
$\mathbf{I}^{(t)}$, and $\mathbf{G}^{(t)}$ the cell state of the previous time
step $\mathbf{C}^{(t-1)}$ is updated as follows

$\mathbf{C}^{(t)}=\mathbf{F}^{(t)}\odot\mathbf{C}^{(t-1)}+\mathbf{I}^{(t)}\odot\mathbf{G}^{(t)},$
(19)

where $\odot$ denotes the element-wise vector (Hadamard) product. The first
term in Eq. (19) removes parts of the previous cell state through the forget
gate and the second term adds new information yielding the new cell state.

Finally, the hidden state $\mathbf{h}^{(t)}$ is updated as follows

$\displaystyle\mathbf{O}^{(t)}$ $\displaystyle=$
$\displaystyle\sigma\left(\mathbf{b}^{O}+\mathbf{W}_{h}^{O}\cdot\mathbf{h}^{(h-1)}+\mathbf{W}_{x}^{O}\cdot\mathbf{x}^{(t)}\right),$
(20) $\displaystyle\mathbf{h}^{(t)}$ $\displaystyle=$
$\displaystyle\mathbf{O}^{(t)}\odot\tanh\left(\mathbf{C}^{(t)}\right),$ (21)

where $\mathbf{O}^{(t)}$ is known as the output gate which uses a sigmoid
activation function for gating and the corresponding bias parameters
$\mathbf{b}^{O}\in\mathbb{R}^{k}$ and weights
$\mathbf{W}_{h}^{O}\in\mathbb{R}^{k\times k}$ and
$\mathbf{W}_{x}^{O}\in\mathbb{R}^{k\times m}$. The new hidden state
$\mathbf{h}^{(t)}$ is then used to compute what to forget, input, and output
by the cell in the next time step. Thus, while both hidden state and cell
state control the memory of the network, the cell state carries information
about the entire sequence and the hidden state encodes the information about
the most recent time step. Similarly to a deep vanilla RNN architecture, in a
deep LSTM NN the hidden state of each LSTM cell is passed to the next layer,
as shown in Fig. 1a.

Since each of the four gates are based on a perceptron the total number of
trainable parameters of the LSTM variant shown in Fig. 1c is $4k(k+m+1)$ and,
since $\mathbf{b}^{F},\mathbf{b}^{I},\mathbf{b}^{O}$, $\mathbf{W}_{h}^{F}$,
$\mathbf{W}_{h}^{I}$, $\mathbf{W}_{h}^{O}$, $\mathbf{W}_{x}^{F}$,
$\mathbf{W}_{x}^{I}$, and $\mathbf{W}_{x}^{O}$ are shared across all LSTM
cells, this is also the total number of trainable parameters of the first LSTM
layer comprised of any number of LSTM cells. The total number of trainable
parameters of a two-layer LSTM architecture is
$4[k_{2}(k_{1}+k_{2}+1)+k_{1}(k_{1}+m+1)]$. Thus, such LSTM models contain
exactly four times more trainable parameters than the vanilla RNN models with
the same number of units (length of the hidden state vector).

####  II.2.5 Gated recurrent unit networks

LSTM has become a popular off-the-shelf architecture that effectively solves
the vanishing gradient problem. However, LSTM is often criticized for its ad
hoc nature. Furthermore, the purpose of its many components is not apparent
and there is no proof that LSTM is even the optimal structure. For example,
Ref. Jozefowicz  _et al._ , 2015 shows that the forget gate is crucial element
of the LSTM architecture, while the output gate is the least important.

Gated recurrent unit (GRU) Cho  _et al._ (2014); Chung  _et al._ (2016, 2015);
Jozefowicz  _et al._ (2015) ANNs are another example of gated RNNs. GRUs were
introduced as a simplified version of LSTM that uses one less gate and, thus,
has fewer trainable parameters. The main difference between LSTM and GRU cells
is that forget and input gates appear in the LSTM architecture are combined
together in one gate in the GRU cell. Additionally, the hidden state and cell
state are combined as well. The information between GRU cells is transferred
via the hidden state vector. Due to the reduction in trainable parameters and
complexity models based on GRU cells tend to converge faster.

There exist several variants of a GRU cell. Fig. 1d illustrates the GRU cell
used in this work. The update equations are discussed below. An input vector
$\mathbf{x}^{(t)}$ is concatenated with the previous cell hidden state
$\mathbf{h}^{(t-1)}$ and passed through the update $Z$ and reset $R$ gates
both using the sigmoid activation function

$\displaystyle\mathbf{Z}^{(t)}$ $\displaystyle=$
$\displaystyle\sigma\left(\mathbf{b}^{Z}+\mathbf{W}_{h}^{Z}\cdot\mathbf{h}^{(t-1)}+\mathbf{W}_{x}^{Z}\cdot\mathbf{x}^{(t)}\right),$
(22) $\displaystyle\mathbf{R}^{(t)}$ $\displaystyle=$
$\displaystyle\sigma\left(\mathbf{b}^{R}+\mathbf{W}_{h}^{R}\cdot\mathbf{h}^{(t-1)}+\mathbf{W}_{x}^{R}\cdot\mathbf{x}^{(t)}\right),$
(23)

where
$\mathbf{b}^{Z(R)}=\mathbf{b}_{W_{h}}^{Z(R)}+\mathbf{b}_{W_{x}}^{Z(R)}\in\mathbb{R}^{k}$
are the biases and $\mathbf{W}_{h}^{Z(R)}\in\mathbb{R}^{k\times k}$,
$\mathbf{W}_{x}^{Z(R)}\in\mathbb{R}^{k\times m}$ are the weights of the update
(reset) gate. The update gate controls the update of the memory state
$\mathbf{h}^{(t)}$. The reset gate controls the influence of the hidden state
$\mathbf{h}^{(t-1)}$ on $\mathbf{G}^{(t)}$ introducing additional nonlinear
effect in the relationship between past and future states. The two gates can
individually “ignore” parts of the hidden state vector. Goodfellow  _et al._
(2016)

Next, the input vector $\mathbf{x}^{(t)}$ is concatenated with the element-
wise product of the hidden state vector $\mathbf{h}^{(t-1)}$ and the output of
the reset gate which, after passing through hyperbolic tangent activation
function, produces the so-called proposed new candidate state

$\mathbf{G}^{(t)}=\tanh\left(\mathbf{b}^{G}+\mathbf{W}_{h}^{G}\cdot\left(\mathbf{R}^{(t)}\odot\mathbf{h}^{(t-1)}\right)+\mathbf{W}_{x}^{G}\cdot\mathbf{x}^{(t)}\right),$
(24)

where
$\mathbf{b}^{G}=\mathbf{b}_{W_{h}}^{G}+\mathbf{b}_{W_{x}}^{G}\in\mathbb{R}^{k}$,
$\mathbf{W}_{h}^{G}\in\mathbb{R}^{k\times k}$, and
$\mathbf{W}_{x}^{G}\in\mathbb{R}^{k\times m}$ are the corresponding bias
vector and weights. The hidden state is updated by combining the previous cell
hidden state $\mathbf{h}^{(t-1)}$ with the output of the update gate as well
as with the product of the output of the update gate and the proposed
candidate gate

$\mathbf{h}^{(t)}=\left(1-\mathbf{Z}^{(t)}\right)\odot\mathbf{h}^{(t-1)}+\mathbf{Z}^{(t)}\odot\mathbf{G}^{(t)}.$
(25)

The total number of trainable parameters of a single GRU cell described above,
as implemented in TensorFlow software library, is $3k(k+m+2)$. Note the use of
the two separate sets of biases $\mathbf{b}_{W_{h}}$ and $\mathbf{b}_{W_{x}}$.
The total number of trainable parameters is independent of the number of GRU
cells due to parameter sharing. Two-layer GRU models contain
$3[k_{2}(k_{2}+k_{1}+2)+k_{1}(k_{1}+m+2)]$ trainable parameters. Thus, given
the same $k_{1}$, $k_{2}$, and $m$, GRU models contain approximately three
times more parameters than the corresponding vanilla RNN models and ca. 1/4
less parameters than LSTM models with the same number of layers.

####  II.2.6 Bidirectional RNNs

Figure 2: Bidirectional recurrent neural networks. Blue rounded rectangles
represent cells and $x_{t}$ is the $t$th element of an input sequence
$\mathbf{x}$. The input sequence $\mathbf{x}$ is processed by the two
disconnected RNNs from beginning to the end and vice versa. Correspondingly,
two hidden state vectors forward $\overrightarrow{\mathbf{h}}^{(t)}$ and
backward $\overleftarrow{\mathbf{h}}^{(t)}$ are updated and passed between the
cells of the same layer.

In all the RNNs, described above the state at time $T$ captures information
from the past $\left(x^{(t)}\right),t=1,\ldots,T-1$ as well as the current
state $x^{(T)}$ to make a prediction. Such RNNs are said to have “causal”
structure Goodfellow  _et al._ (2016) and are called unidirectional RNNs. RNN
architectures that combine an RNN that moves forward through time from $t=1$
to $t=T$ with another RNN that moves backward through time from $t=T$ to $t=1$
are called bidirectional RNNs. Graves  _et al._ (2013); Goodfellow  _et al._
(2016) This is realized by duplicating each recurrent layer in the network.
The two resulting layers have separate forward $\overrightarrow{\mathbf{h}}$
and backward $\overleftarrow{\mathbf{h}}$ hidden state vectors. Forward and
backward layers are not connected to each other. An input time sequence is
provided in the chronological order to the first RNN layer and is fed in the
reversed chronological order to the second RNN layer. Applying RNNs twice
increases the amount of input information available to the network and leads
to better capturing long-term dependencies and, thus, improves the accuracy of
the model. Baldi  _et al._ (1999)

Bidirectional recurrent NN can be built from the vanilla RNN, LSTM, and GRU
cells resulting in BRNN, BLSTM, Schuster and Paliwal (1997) BGRU models, Lynn
_et al._ (2019) respectively. In each case the update equations are similar to
Eqs. (15)-(25) but applied separately to forward and backward hidden state
vectors. For example, for BRNN with tanh activation function,
$\overrightarrow{\mathbf{h}}$ and $\overleftarrow{\mathbf{h}}$, assuming they
have the identical size, $\overrightarrow{\mathbf{h}}$,
$\overleftarrow{\mathbf{h}}\in\mathbb{R}^{k}$, are updated as follows

$\displaystyle\overrightarrow{\mathbf{h}}^{(t)}$ $\displaystyle=$
$\displaystyle\tanh\left(\mathbf{b}_{\overrightarrow{h}}+\mathbf{W}_{\overrightarrow{h}}\cdot\overrightarrow{\mathbf{h}}^{(t-1)}+\mathbf{W}_{x\overrightarrow{h}}\cdot\mathbf{x}^{(t)}\right),$
(26) $\displaystyle\overleftarrow{\mathbf{h}}^{(t)}$ $\displaystyle=$
$\displaystyle\tanh\left(\mathbf{b}_{\overleftarrow{h}}+\mathbf{W}_{\overleftarrow{h}}\cdot\overleftarrow{\mathbf{h}}^{(t+1)}+\mathbf{W}_{x\overleftarrow{h}}\cdot\mathbf{x}^{(t)}\right),$
(27)

where
$\mathbf{b}_{\overrightarrow{h}},\mathbf{b}_{\overleftarrow{h}}\in\mathbb{R}^{k}$
are the bias parameters for the forward and backward hidden states,
respectively;
$\mathbf{W}_{\overrightarrow{h}},\mathbf{W}_{\overleftarrow{h}}\in\mathbb{R}^{k\times
k}$ and
$\mathbf{W}_{x\overrightarrow{h}},\mathbf{W}_{x\overleftarrow{h}}\in\mathbb{R}^{k\times
m}$ are the coefficient matrices. In deep BRNN architectures, forward and
backward hidden states of a previous layer are combined

$\mathbf{y}^{(t)}=\mathcal{F}\left(\overrightarrow{\mathbf{h}}^{(t)},\overleftarrow{\mathbf{h}}^{(t)}\right),$
(28)

and passed to the following layer. In Eq. (28) $\mathcal{F}$ is a function
that combines the two hidden state vectors. It can be a concatenating function
(this work), element-wise addition, multiplication, or averaging. Chollet  _et
al._ (2015); Abadi  _et al._ (2015) Thus every hidden RNN layer receives and
input from both forward and backward layers at the preceding layer.

The total number of trainable parameters of the first layer of BRNN is exactly
twice the number of trainable parameters in the first layer of a
unidirectional RNN with the same length of the hidden state vector and same
type of RNN cell. In a deep bidirectional recurrent NN architecture an $i$th
hidden layer with $m$ RNN cells and the length of hidden vectors $k_{i}$
outputs to the next RNN layer a tensor $\mathbf{Z}\in\mathbb{R}^{2k_{i}\times
m}$, assuming $\mathcal{F}$ is a concatenating function. The next $(i+1)$
layer uses $2k_{i+1}(2k_{i}+k_{i+1}+1)$ trainable parameters, where $k_{i+1}$
is the length of the hidden state vector of $(i+1)$ layer. Therefore, deep
bidirectional recurrent NNs, in general, use more than twice the number of
trainable parameters than the corresponding unidirectional RNN with the same
number of cells and length of the hidden state vector.

####  II.2.7 Convolutional RNNs

Convolutional Recurrent Neural networks are deep ANNs that combine
convolutional layers with recurrent layers. In this work, we build two-layer
convolutional recurrent NN models by adding one recurrent layer (vanilla RNN,
GRU, LSTM) to 1D CNN layer resulting in three models denoted as CRNN, CGRU,
and CLSTM. The update equations of each of these types of neural networks are
exactly those of the corresponding individual layers described above. For
example, in the CRNN architecture, the first layer (1D CNN) processes input
sequence of length $T$ through convolution, or rather, cross-correlation,
operation as given by Eqs. (10) and (11) and generates output tensor
$\mathbf{Z}^{(1)}$ formed by the output of $K_{1}$ kernels (filters) each of
size $M_{1}=\lfloor(T-k_{1}+2p_{1})/s_{1}\rfloor+1$, where $s_{1}$ is the
stride and $p_{1}$ is the zero padding. The output tensor $\mathbf{Z}^{(1)}$
is then passed to an RNN layer which is comprised of $M_{1}$ vanilla RNN cells
with the hidden state vector of (user-specified) length $k_{2}$. Each RNN cell
receives an input $x^{(t)}\equiv Z^{(1)}_{t,1:K_{1}},i\in 1,\ldots,M_{1}$ from
the preceding 1D CNN layer and processes it according to the corresponding
update equations.

In this work the recurrent layer of convolutional recurrrent NN models is
configured to output the whole sequence of hidden states from all RNN cells,
$\mathbf{Z}^{(2)}\in\mathbb{R}^{M_{1}\times k_{2}}$ which is then passed to
fully-connected layers. The total number of trainable parameters of the two-
layer architecture described above is $K_{1}(k+1)+N_{\mathrm{RNN}}$, where
$k_{1}$ is the size of the kernel of 1D CNN layer and $N_{\mathrm{RNN}}$ the
number of trainable parameters of an RNN layer.

####  II.2.8 Convolutional bidirectional RNNs

It is possible to fuse a 1D CNN layer and a bidirectional RNN layer. Resulting
architectures are called Convolutional Bidirectional Recurrent Neural
Networks. Such models are expected to combine the benefits of convolutional
recurrent NNs with the improved description of long-term dependencies
pertinent to bidirectional RNNs. In this work, three convolutional
bidirectional NN models are build by combining one 1D CNN layer with each of
the three recurrent NN architectures discussed above: the vanilla RNN, GRU,
and LSTM. The resulting models are denoted as CBRNN, CBGRU, and CBLSTM.
Similarly to convolutional recurrent NNs the update equations of convolutional
bidirectional NNs can be deduced from the corresponding update equations of 1D
CNN and bidirectional recurrent NNs given above. Specifically for CBRNN, an
input sequence is first passed through an 1D CNN layer where it is transformed
according to Eqs. (10) and (11) into an output tensor $\mathbf{Z}^{(1)}$ which
is formed by the output of $K_{1}$ kernels (filters) each of size
$M_{1}=\lfloor(T-k_{1}+2p_{1})/s_{1}\rfloor+1$, where $s_{1}$ is the stride
and $p_{1}$ is the zero padding. Then, $\mathbf{Z}^{(1)}$ is fed separately
into forward and backward vanilla RNN layers where the corresponding forward
and backward hidden states are updated as shown in Eqs. (26) and (27). In the
two-layer architecture described above, the total number of trainable
parameters is simply the sum of the total number of trainable parameters of
the 1D CNN layer and the total number of trainable parameters of a
bidirectional recurrent layer.

####  II.2.9 Kernel ridge regression

In kernel ridge regression (KRR) Stulp and Sigaud (2015); Dral (2019a); Hastie
_et al._ (2009) the approximating function $f(\mathbf{x})$ for a vector of
input values $\mathbf{x}$ is defined as

$f(\mathbf{x})=\sum_{i=1}^{N}\alpha_{i}k\left(\mathbf{x},\mathbf{x}_{i}\right),$
(29)

where $N$ is the number of training points and $\bm{\alpha}=\\{\alpha_{i}\\}$
is a vector of regression coefficients. The covariance function
$k\left(\mathbf{x},\mathbf{x}_{i}\right)$, commonly referred to as a kernel
function or simply kernel, can be understood as a similarity measure between
two vectors $\mathbf{x}$ and $\mathbf{x}_{i}$ from the input space. The kernel
performs an implicit mapping to a higher-dimensional feature space.

Because input time-sequences employed in this work have the same length we
focus on standard kernels. Gneiting  _et al._ (2010) One of the most common
kernel functions is the Matérn kernel Rasmussen and Williams (2005); Gneiting
_et al._ (2010), which in the MLatom software package Dral (2019b); Dral  _et
al._ (2021, 2022) used in this work is defined as: Dral  _et al._ (2021)

$\displaystyle k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}}{\sigma}\right)\sum_{k=0}^{n}\frac{(n+k)!}{(2n)!}$
(30) $\displaystyle\times$
$\displaystyle\binom{n}{k}\exp\left(\frac{2\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}}{\sigma}\right)^{n-k},$

where $\sigma$ is a positive hyperparameter which defines the characteristic
length scale of the covariance function, $n$ is a non-negative integer, and
$\|\ldots\|_{2}$ is the Euclidian distance which is taken to be the $L^{2}$
norm. For $n=0$ the Matern covariance function reduces to exponential kernel
function Rasmussen and Williams (2005); Gneiting  _et al._ (2010)

$k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)=\exp\left(-\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}}{\sigma}\right).$
(31)

Another popular choice of a covariance function is the squared exponential
(Gaussian) kernel function Rasmussen and Williams (2005)

$k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)=\exp\left(-\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}_{2}}{2\sigma^{2}}\right).$
(32)

Because quantum dynamics often resembles periodically-decaying time-series, we
also test a decaying periodic kernel (our adaptation based on Refs. 153; 157):

$k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)=\exp\Bigg{(}-\frac{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}_{2}}{2\sigma^{2}}-\frac{2}{\sigma_{p}^{2}}\sin^{2}{\left(\frac{\pi}{p}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|_{2}\right)}\Bigg{)},$
(33)

where $p$ is the period and $\sigma_{p}$ is a length scale for the periodic
term (both are hyperparameters).

Given the kernel function, the regression coefficients $\bm{\alpha}$ are found
by minimizing a squared error loss function

$\min_{\alpha}\sum_{i=1}^{N}\left(f\left(\mathbf{x}_{i}\right)-y_{i}\right)^{2}+\lambda\bm{\alpha}^{T}\mathbf{K}\bm{\alpha},$
(34)

where $\mathbf{y}=\\{y_{i}\\}$ is the target output vector,
$\mathbf{K}\in\mathbb{R}^{N\times N}$ is the kernel matrix with elements
$K_{ij}=k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$ and $\lambda$ denotes a
non-negative regularization hyperparameter. In Eq. (34), the second term is
usually added to prevent KRR model from assigning large weight to a single
point. The optimization of parameters (regression coefficients $\bm{\alpha}$)
amounts to solving a system of linear equations

$\left(\mathbf{K}+\lambda\mathbf{I}\right)=\mathbf{y},$ (35)

for which analytical solution is known. Here $\mathbf{I}$ is the identity
matrix. The computational scaling of solving this system is
$\mathcal{O}(N^{3})$. Dral  _et al._ (2021); Rasmussen and Williams (2005)

KRR is a kernelized version of ridge regression and when linear kernel
function

$k\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)=\mathbf{x}_{i}^{T}\mathbf{x}_{j}$
(36)

is used, KRR becomes equivalent to ridge regression, i.e., the approximating
function $f(\mathbf{x})$ is simply a multiple linear regression with
regression coefficients $\bm{\beta}$ shrunk (see the second term in Eq. (34))
using regularization:

$\displaystyle f(\mathbf{x})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\alpha_{i}\mathbf{x}_{i}^{T}\mathbf{x}=\sum_{i=1}^{N}\alpha_{i}\sum_{s=1}^{T}x_{is}x_{s}$
(37) $\displaystyle=$
$\displaystyle\sum_{s=1}^{T}\left(x_{s}\sum_{i=1}^{N}\alpha_{i}x_{is}\right)=\sum_{s=1}^{T}\beta_{s}x_{s}.$

As shown in above equation, regression coefficients $\bm{\beta}$ can be
conveniently derived from $\bm{\alpha}$ coefficients and training input data
and, in fact, are printed out by MLatom.

##  III Computational details

###  III.1 Data sets for training, validation, and testing

The data set used in this work is the same as used in Ref. Ullah and Dral,
2021 and can be accessed at https://doi.org/10.6084/m9.figshare.15134649. It
was generated as detailed below. Ullah and Dral (2021) Firstly, HEOM
calculations for all combinations of the following system and bath parameters:
$\epsilon/\Delta=\\{0,1\\}$,
$\lambda/\Delta=\\{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0\\}$,
$\omega_{c}/\Delta=\\{1,2,3,4,5,6,7,8,9,10\\}$, and
$\beta\Delta=\\{0.1,0.25,0.5,0.75,1\\}$, were performed with QuTiP software
package. Johansson  _et al._ (2012) In our calculations we set $\Delta=1.0$.
The total propagation time was $t_{\mathrm{max}}\Delta=20$ and the HEOM
integration time-step was set to $t\Delta=0.05$. In total, 1,000 HEOM
calculations, 500 for symmetric ($\epsilon/\Delta=0$) and 500 for asymmetric
($\epsilon/\Delta=1$) spin-boson Hamiltonian, were performed. Time-evolved
reduced density matrices (RDM) are saved every $dt\Delta=0.1$. Secondly,
$\langle\hat{\sigma}_{z}(t)\rangle$ are calculated from RDMs and processed
into shorter sequences of length $T$ by window slicing. Herrera Rodríguez and
Kananenka (2021); Ullah and Dral (2021); Lin  _et al._ (2021) Namely, for a
time series $\mathbf{x}=\left(x^{(1)},\ldots,x^{(L)}\right)$, where
$\langle\hat{\sigma}_{z}(t)\rangle$ is denoted by $x^{(t)}$ for compactness, a
slice is a subset of the original time series defined as
$\mathbf{s}_{i:j}=\left(x^{(i)},\ldots,x^{(j)}\right),1\leq i\leq j\leq P$.
For a given time series $\mathbf{x}$ of length $L$, and the length of the
slice is $P$, a set of $L-P+1$ sliced time series
$\\{\mathbf{s}_{1:P},\mathbf{s}_{2:P+1},\ldots,\mathbf{s}_{L-P+1:L}\\}$ is
generated. Finally, the total data set
$\mathcal{D}=\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{N}$ containing time series
$\mathbf{x}_{i}$ and their corresponding labels $y_{i}$ is obtained by setting
$1,\ldots,T$ elements of each slice, with $T=P-1$, to an input time-series
$\mathbf{x}_{i}$ and the last ($P$th) element of each slice to the associated
label $y_{i}$.

In general, the size of the window $P-1$, or equivalently $T$, should be
treated as a hyperparameter but, following previous work, Ullah and Dral
(2021) we set to $T=0.2L$. The window slicing is applied to all 1,000 RDMs
obtained in HEOM calculations with different system and system-bath
parameters. For each set of the Hamiltonian parameters the initially
calculated set of time-evolved $\langle\hat{\sigma}_{z}(t)\rangle$ with
$L=t_{\mathrm{max}}/dt=200$ generates 160 data points for the data set with
$T=41$ (including $t\Delta=0$ point).

From the raw HEOM data set of 1,000 trajectories, 100 randomly chosen
trajectories are taken as the hold-out test set, which is used for testing and
generating the results presented in Sec. IV. The remaining set of 900
trajectories are transformed into 144,000 trajectories by window slicing
described above. In total 72,000 short-time
$\langle\hat{\sigma}_{z}(t)\rangle$ trajectories for symmetric and 72,000 for
asymmetric spin-boson models were generated. Each trajectory has a time-length
of $t\Delta+dt\Delta=4.1$. This data set of supervised trajectories is the
training set which is randomly partitioned into two subsets: a sub-training
set, which contains 80% of the data and a validation set with 20% of the data.
ML models are (initially) trained on the sub-training set and the validation
set is used for monitoring the performance of the models (mainly, to prevent
overfitting and optimize hyperparameters). The final NN models tested on a
hold-out test set are not trained on the entire training set and are only
trained on the sub-training set to prevent overfitting. KRR models are,
however, trained on the entire training set. The performance of KRR models
trained only on the sub-training set is similar to the KRR models trained on
the entire training set (Appendix A). These are typical ways of training NN
and KRR models, respectively, which are different due to differences in the
formalism and training procedures of these types of models. Following previous
similar works Herrera Rodríguez and Kananenka (2021); Ullah and Dral (2021)
the input data is not normalized.

###  III.2 Artificial neural network models

####  III.2.1 Details of the models

In this work fourteen deep ANN models are built and tested. In general, each
ANN model is comprised by the total of two convolutional and/or recurrent
layers followed by one fully-connected layer, and one output layer. The
exception is the FFNN model which comprised of two fully-connected layers
followed by the output layer. The fully-connected layer in each model, except
for FFNN, has 256 neurons and the rectified linear function defined as (ReLU)
$f(z)=\max(0,z)$ is used as the activation function. Fixing the properties of
fully-connected layer allows to compare the performance of recurrent and
convolutional layers. The output layer contains one neuron with the linear
activation function $f(z)=z$. The details of all ANN models studied in the
present Article is summarized below.

  * •

1D CNN model contains two 1D CNN layers followed by a MaxPooling, one fully-
connected and output layers. ReLU is used as the activation function in each
1D CNN and fully-connected layers. For the MaxPooling layer a pool size
$k_{p}=2$ (see Eq. (13)) is used. The stride of $s=1$ and zero padding are
used in both 1D CNN and MaxPooling layers. The number of filters and the
filter sizes of each layer are optimized using Particle Swarm Optimization
algorithm as described in Sec. III.2.2.

  * •

FFNN model contains two hidden fully-connected layers followed by the output
layer. The ReLU activation function is used.

  * •

Three recurrent NN models comprised of the two recurrent layers of the same
type: the vanilla RNN, LSTM, and GRU. The whole sequence of hidden state
vectors from all cells is passed from the first to the second recurrent layer
as well as from the second recurrent layer to the fully-connected layer.

  * •

Three bidirectional recurrent NN models comprised of the two bidirectional
recurrent layers of the same type: the vanilla BRNN (denoted simply as BRNN
hereafter), BLSTM, and BGRU. The whole sequence of hidden state vectors from
all cells is passed from the first bidirectional recurrent layer to the second
bidirectional recurrent layer as well as from the second bidirectional
recurrent layer to the fully-connected layer.

  * •

Three convolutional recurrent NN models comprised of one 1D CNN layer and one
recurrent layer of each type: the vanilla RNN, LSTM, and GRU. The resulting
models are denoted as CRNN, CLSTM, and CGRU. The recurrent layer returns the
whole sequence of hidden state vectors from all cells to the fully-connected
layer. The stride of $s=1$, zero padding, and the ReLU activation function are
used in the 1D CNN layer. The number of filters and the filter size is the
same as in the first layer of the 1D CNN model.

  * •

Three convolutional bidirectioinal recurrent NN models comprised of one 1D CNN
layer and one bidirectional recurrent layer of each type: the vanilla RNN,
LSTM, and GRU. The resulting models are denoted as CBRNN, CBLSTM, and CBGRU.
The stride of $s=1$, zero padding, and the ReLU activation function is used in
the 1D CNN layer. The number of filters and the filter size is the same as in
the first layer of the 1D CNN model. The bidirectional recurrent layer returns
the whole sequence of the forward-backward hidden state vectors from all cells
to the fully-connected layer.

####  III.2.2  Hyperparameter optimization

The goal of this work is to compare the performance of ML models including ANN
models with different architectures as described above. Depending on the
overall objective, one can envision several approaches for comparing
performances of different ANNs. For example, the hyperparameters of each of
ANN model can be adjusted using grid search, random search, or other
optimization techniques, including evolutionary algorithms, to achieve the
best performance of each ANN model on the given data set. However, given the
fundamentally different types of ANN models considered in the present Article,
this approach may result in several models with approximately the same
prediction accuracy but requiring drastically different resources such as CPU
(or GPU) time and memory. In such case, it might be reasonable to use ANN
models that are not the most accurate, but provide an acceptable accuracy with
less training and prediction times. The notion of the optimal running time,
however, strongly depends on the problem under study. One can also set the
desired accuracy level and systematically adjust the hyperparameters of each
model to achieve the desired accuracy and then analyze the resulting ANN
models. This approach, however, requires specifying the fixed accuracy level
which might not be achievable for some ML models and relies on a priori
knowledge of the dynamics of the system of interest.

The following approach is adopted in the present study. We compare the
performance of ANN models with approximately the same number of trainable
parameters. Hyperparameters of each ANN model described above are scanned
using a grid search approach and the models with the pre-determined number of
trainable parameters are selected and their performance is reported in Sec.
IV. The number of trainable parameters is chosen as follows. In Refs. 32 and
159 we illustrated that deep 1D CNN models can approximate long-time dynamics
of a molecular dimer system (spin-boson-like model) and the
Fenna–Matthews–Olson photosynthetic complex accurately. Therefore, the 1D CNN
model from Ref. 32 is taken as a base model. The hyperparameters of this
model, associated with convolutional layers, specifically, the number and the
size of the kernels of each layer, are optimized using Particle Swarm
Optimization algorithm (PSO). All other hyperparameters such as the number of
neurons of a fully-connected layer and the activation functions of all neurons
are fixed to 256 and ReLU correspondingly. PSO calculation is performed with
three particles and for 50 steps. The data set described above is used for the
hyperparameter optimization. At each step of the PSO calculation, three 1D CNN
models, one for each particle, are trained and validated using the above
mentioned data set. During training, the deviation between predicted
$\hat{y}_{i}$ and reference values $y_{i}$ is minimized. In this work, we use
the mean squared error (MSE) as the loss function in the minimization

$\mathrm{MSE}=\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}\left(y_{i}-\hat{y}_{i}\right)^{2}.$
(38)

Adaptive moment estimation (Adam) algorithm Kingma and Ba (2017) is used with
the initial learning rate set to 1.0$\cdot$10-4. The initial values of the
weights are randomly sampled using Xavier initialization. Glorot and Bengio
(2010) The biases are initialized to zero. Keras Chollet  _et al._ (2015)
software package with the TensorFlow Abadi  _et al._ (2015) backend was
employed for calculation. The batch size is set to $N_{b}=64$ and the training
for each model is performed for 30 epochs. As will be shown later this is
sufficient for the MSE to drop below 10-5–10-6.

For each trained 1D CNN model, the MSE for the validation set is calculated
and used as the fitness value in the PSO. Then the hyperparameters of 1D CNN
models are adjusted based on the algorithmic details of PSO which can be found
in Appendix B. After 50 steps of the PSO, the 1D CNN models stop improving and
the PSO calculation is terminated. The optimized number of kernels of the
final 1D CNN model is 235 and 125 with the kernel sizes 16 and 7 of the first
and second layer, respectively. The 1D CNN model described above contains a
total of 530,258 trainable parameters. This number is taken as reference and
the hyperaprameters of all other 13 ANN models are selected to generate models
with approximately the same number of parameters. Note, however, that
according to the formulas given in Sec. I, it is not possible to set
hyperparameters in all models to achieve exactly the desired number of
trainable parameters. Therefore, some variation in the number of trainable
parameters across the reported ANN models is to be expected.

Even though the number of trainable parameters is straightforwardly connected
to the hyperparameters of each ANN architecture, a grid search is performed
over hyperparameters with the goal to understand the sensitivity of the model
performance to the small variations in the hyperparameters. In all models
containing (bidirectional) recurrent layers the number of units is scanned
from 5 to 100 with the step of 5. Additionally, in all ANN models containing
1D CNN layers and (bidirectional) recurrent layers, the number of kernels and
kernel sizes of the 1D CNN layer is taken to be same as in 1D CNN model while
the number of units in the recurrent layer is scanned from 5 to 100 in steps
of 5. All the models are trained, validated, and tested using the data set for
the symmetric spin-boson model. keras software package and Adam algorithm is
used for all ANN calculations. The learning rate is fixed to 1.0$\cdot$10-4,
the batch size is set to 128, and the training for each model is performed for
30 epochs.

Usually more than one set of hyperparameters generates ANN models with the
number of trainable parameters close to the target number. In such cases,
models with the number of trainable parameters within $\pm$5% of the reference
number are considered. Among selected models, the models that perform
significantly different than the average model are discarded. The number of
outliers in each case is found to be small. We attribute outliers to finite
training time and nonuniform decay of the MSE during the training, which in
case of RNNs is a manifestation of the vanishing gradient problem. We will
return to this issue in Sec. IV. The model with the closest number of
trainable parameters to the target number is saved and its performance is
reported.

####  III.2.3 Training and validation curves

To illustrate the learning process of ANN models in Fig. 3 we plot the
training and validation curves for several such models. Such curves are a
widely used tool to examine the performance of supervised learning algorithms.
In general the learning is found to be stable in each case. The apparent
noisiness in the data should be attributed to the logarithmic scale on the
vertical axis. Several conclusions can be drawn by examining the curves shown
in Fig. 3. Firstly, none of the models overfits which is illustrated in lower
and in, general, decaying validation MAE (Fig. 3 right panel) compared to the
training MAEs (Fig. 3 left panel). Secondly, the learning process is fast even
with the learning rate of 1$\cdot 10^{-4}$. It should be noted that decreasing
learning rate to 1$\cdot 10^{-5}$ does not result in noticeable improvement in
accuracy for all the ANN models studied in this work but the increase of the
learning rate to 1$\cdot 10^{-3}$ usually results in increasing MAE.
Furthermore, FFNN model reaches the lowest MAE rapidly in less than 10
training epochs while the MAE of other models e.g., 1D CNN, B(LSTM,GRU)
continues to decay over all shown 30 training epochs. The same behavior is
manifested in the validation MAEs. This suggests that 30 training epochs used
in this work is justified for 1D CNN and other convolutional and recurrent
models but seems unnecessarily too many for the FFNN model. This observation
further illustrates the efficiency of simple FFNN models. We note that in the
case of the BRNN model the validation MAE exhibits large-magnitude changes
over the course of training. This can be attributed to the vanishing gradient
problem of the vanilla RNNs.

Figure 3: Learning curves for selected ANN models: training (left) and
validation (right). Note the logarithmic scale of the vertical axis.

###  III.3 KRR models

In addition to 14 ANN models, 8 KRR models are investigated in the present
Article. Each KRR model differs by the choice of the kernel. The following
kernels are studied: the linear kernel denoted as KRR-L (Eq. (36)), Gaussian
kernel (Eq. (32), KRR-G), exponential kernel Eq. (31) (KRR-E), Matérn kernel,
Eq. (30) with $n=1,2,3,4$ (KRR-M1, KRR-M2, KRR-M3, KRR-M4), and a periodic-
decaying kernel, Eq. (33) (KRR-DP). It should be noted that unlike recurrent
neural networks, KRR (and FFNN) models studied in this work require a fixed-
size input. To enable the comparison between ANN and KRR models studied in
this work, we focused on the input $\langle\hat{\sigma}_{z}(t)\rangle$
trajectories of the same length. We note that extending RNN models built in
the present Article to variable size input is straightforward and will be
discussed elsewhere.

In KRR models, optimization amounts to finding regression coefficients
$\bm{\alpha}$ as shown in Eq. (35). Therefore, one can interpret each
regression coefficient $\alpha_{i}$ as a trainable parameter. The total number
of such coefficients is the same as the number of elements of the training
set, which in this work is 72,000 in each case of symmetric and asymmetric
spin-boson models. This number of trainable parameters of KRR model is far
fewer than the target number of trainable parameters for ANN models. The
discrepancy between the number of trainable parameters which, by construction,
is set to the size of the training set, yet it is variable in ANN models. This
further makes the faithful comparison of fundamentally different ML approaches
studied in this work non-trivial. As will be shown in the following, the fewer
number of trainable parameters of KRR does not adversely impact the
performance of KRR models.

KRR approaches have analytical solution to optimal parameters (regression
coefficients $\bm{\alpha}$), but they still require adjusting a (small) number
of hyperparameters. All KRR models used here contain the regularization
hyperparameter $\lambda$ which is the only hyperparameter in KRR-L. Other KRR
models have additional hyperparameters in their kernel functions, i.e., all
nonlinear models have length-scale hyperparameter $\sigma$, KRR-M$n$ also have
the integer hyperparameter $n$ (not optimized here), KRR-DP has two more
hyperparameters compared to KRR-G: $p$ and $\sigma_{p}$. Due to the small
number of hyperparameters, they are optimized on a logarithmic grid as
described elsewhereDral (2019b) for all models except for KRR-DP which has
three hyperparameters. Hyperparameters in KRR-DP are optimized using the tree-
structured Parzen estimator Bergstra  _et al._ (2011) via interface to the
hyperopt package Bergstra  _et al._ (2015). All hyperparameters of KRR models
are optimized based on the data set for symmetric spin-boson model using grid-
search method for KRR.

The hyperparameters of each ML model are optimized only for the data set for
symmetric spin-boson model as already indicated above. The optimized models
are then re-trained on the asymmetric spin-boson data set keeping the
hyperparameters unchanged. MLatom software package Dral (2019b); Dral  _et
al._ (2021) is used for all KRR calculations.

##  IV Results and discussion

###  IV.1 Symmetric spin-boson model

The details and performance of all ANN and KRR models are summarized in Table
1 and Table 2, respectively. As seen from Table 1, the two best performing ANN
models on the symmetric spin-boson data set, are the convolutional LSTM
(CLSTM) and convolutional bidirectional LSTM (CBLSTM) while the two least
accurate ANN models are unidirectional RNN and bidirectional vanilla RNN
models. Worse performance of simple vanilla RNN layers compared to LSTM and
GRU layers is the expected result. It has been observed in other applications
that in contrast to LSTM and GRU models, the performance of vanilla RNN models
degrades with the increasing length of input sequence. Stérin  _et al._
(2017); Gupta and Saini (2020) We note that the weak performance of RNN models
is despite the increased number of units or length of hidden state vector
compared to other recurrent NN models.

Table 1: Hyperparameters, mean absolute prediction errors, total number of trainable parameters, training, and average single step prediction times of all ANN models studied in this work. For each recurrent layer the number of units is shown. The number of kernels X and kernel sizes Y for each convolutional layer is shown in parenthesis as (X,Y). | Trainable | Layers | Mean absolute error | Time [s]
---|---|---|---|---
Model | Parameters | Layer 1 | Layer 2 | Symmetric | Asymmetric | Training | Prediction
1D CNN | 530,258 | (235,16) | (125,7) | 1.55$\cdot 10^{-3}$ | 4.84$\cdot 10^{-2}$ | 465 | 3.8
FFNN | 520,045 | 754 | 646 | 1.32$\cdot 10^{-3}$ | 3.70$\cdot 10^{-2}$ | 82 | 3.3
Recurrent Neural Networks |
LSTM | 528,577 | 15 | 49 | 1.58$\cdot 10^{-3}$ | 2.35$\cdot 10^{-2}$ | 623 | 6.1
GRU | 553,453 | 60 | 50 | 2.05$\cdot 10^{-3}$ | 2.57$\cdot 10^{-2}$ | 668 | 4.4
RNN | 535,468 | 65 | 50 | 3.00$\cdot 10^{-3}$ | 6.17$\cdot 10^{-2}$ | 302 | 4.2
Convolutional Recurrent Neural Networks |
CLSTM | 501,965 | (28,16) | 71 | 1.17$\cdot 10^{-3}$ | 2.50$\cdot 10^{-2}$ | 279 | 4.1
CGRU | 515,806 | (55,16) | 73 | 1.38$\cdot 10^{-3}$ | 2.14$\cdot 10^{-2}$ | 294 | 4.7
CRNN | 513,673 | (243,16) | 73 | 1.46$\cdot 10^{-3}$ | 3.61$\cdot 10^{-2}$ | 197 | 5.4
Convolutional Bidirectional Recurrent Neural Networks |
CBLSTM | 568,022 | (109,16) | 39 | 1.17$\cdot 10^{-3}$ | 2.84$\cdot 10^{-2}$ | 333 | 3.8
CBGRU | 514,860 | (55,16) | 37 | 1.54$\cdot 10^{-3}$ | 3.68$\cdot 10^{-2}$ | 325 | 5.3
CBRNN | 508,842 | (297,16) | 36 | 2.50$\cdot 10^{-3}$ | 3.58$\cdot 10^{-2}$ | 256 | 3.9
Bidirectional Recurrent Neural Networks |
BLSTM | 511,809 | 6 | 24 | 2.12$\cdot 10^{-3}$ | 2.56$\cdot 10^{-2}$ | 635 | 5.2
BGRU | 534,991 | 14 | 25 | 2.28$\cdot 10^{-3}$ | 2.48$\cdot 10^{-2}$ | 1109 | 6.7
BRNN | 511,959 | 37 | 24 | 6.95$\cdot 10^{-3}$ | 4.27$\cdot 10^{-1}$ | 396 | 5.4

A better performance of convolutional recurrent NN models is an illustration
of the power of sub-sampling of the initial input data which is then more
efficiently learned by the LSTM layer. Comparing the unidirectional or
bidirectional recurrent NN models whether taken separately or as the second
layer on top of 1D CNN layer confirms the trend. The best performing model is
an (C,B,CB)LSTM model, followed by the (C,B,CB)GRU, and (C,B,CB)RNN. It is
also worth noting that recurrent NN models based on bidirectional layers do
not outperform their unidirectional counterparts. This is observed for both
unidirectional recurrent models compared to bidirectional recurrent models as
well as when convolutional unidirectional recurrent models are compared to the
corresponding convolutional bidirectional recurrent models.

Table 2: Mean absolute prediction errors (MAEs), training, and average single step prediction times of all KRR models used in this work. KRR-L, KRR-G, KRR-DP, KRR-E, and KRR-M ($n=1,2,3,4$) denote kernel ridge regression models with linear kernel, Gaussian kernel, decaying-periodic kernel, exponential kernel, and Matern kernel with $n=1,2,3,4$, respectively. | Trainable | Mean absolute error | Time [s]
---|---|---|---
Model | Parameters | Symmetric | Asymmetric | Training | Prediction
KRR-L | 72,000 | 1.2$\cdot 10^{-2}$ | 6.5$\cdot 10^{-2}$ | 196 | 1.6
KRR-G | 72,000 | 4.7$\cdot 10^{-4}$ | 1.2$\cdot 10^{-3}$ | 220 | 1.6
KRR-DP | 72,000 | 4.3$\cdot 10^{-4}$ | 2.0$\cdot 10^{-3}$ | 257 | 1.4
KRR-E | 72,000 | 2.1$\cdot 10^{-3}$ | 3.3$\cdot 10^{-3}$ | 222 | 1.6
KRR-M1 | 72,000 | 2.4$\cdot 10^{-4}$ | 1.3$\cdot 10^{-3}$ | 260 | 1.6
KRR-M2 | 72,000 | 2.2$\cdot 10^{-4}$ | 2.7$\cdot 10^{-3}$ | 259 | 1.7
KRR-M3 | 72,000 | 2.0$\cdot 10^{-4}$ | 2.3$\cdot 10^{-3}$ | 279 | 1.7
KRR-M4 | 72,000 | 2.3$\cdot 10^{-4}$ | 2.1$\cdot 10^{-3}$ | 273 | 1.7

Focusing on KRR models, whose details and performance is illustrated in Table
2, we note that these models generally outperform ANN models with the
exception of a KRR model with linear kernel. Furthermore, the KRR-L is the
worst performing of all ML models studied in the present work. This is,
perhaps, not surprising because it is expected that quantum dynamics of such a
complex system as the spin-boson model is highly non-trivial and cannot be
captured by simple linear regression.

Figure 4: Expectation values $\langle\hat{\sigma}_{z}(t)\rangle$ for symmetric
spin-boson model as a function of time. Results predicted by an indicated ML
model are compared to the HEOM results (black triangles) for the following
parameters: (a) $\lambda=0.2,\omega_{c}=8.0,\beta=1.0$; (b)
$\lambda=0.4,\omega_{c}=10.0,\beta=1.0$; (c)
$\lambda=0.2,\omega_{c}=10.0,\beta=0.25$; (d)
$\lambda=0.1,\omega_{c}=4.0,\beta=0.1$; (e)
$\lambda=0.8,\omega_{c}=3.0,\beta=1.0$; (f)
$\lambda=1.0,\omega_{c}=2.0,\beta=0.1$. All parameters are in the units of
$\Delta$.

The differences in accuracy between ML models noted above is, however, not
significant in the case of symmetric spin-boson model. This is illustrated in
Fig. 4 where the exact population difference
$\langle\hat{\sigma}_{z}(t)\rangle$ is compared to the ML-predicted
$\langle\hat{\sigma}_{z}(t)\rangle$. We stress that only short
$\langle\hat{\sigma}_{z}(t)\rangle$ trajectory of $t\Delta=4$ is used as an
input. The rest of the dynamics is predicted recursively as done in Refs. 32
and 33. Fig. 4 compares the performance of the 1D CNN model, whose number of
trainable parameters is used as a reference for other ANN models, to the best
and the worst performing ANN and KRR models. One notices that prediction
accuracy of the KRR-L model degrades slowly over time but still remains
acceptable. This can only be seen in Fig. 4a and b where the chosen system and
bath parameters generate the oscillatory dynamics of RDM. In the cases of
incoherent relaxation dynamics KRR-L results are indistinguishable from the
exact HEOM dynamics.

We stress that even the worst performing ML model, KRR-L, as illustrated in
Fig. 4 still provides the acceptable accuracy. Therefore, we conclude that all
ML models benchmarked in the present article can be trained to provide a
sufficient long-time prediction accuracy. The suspected drawback of the
recursive propagation approach proposed in Refs. 32 and 33 and exploited in
this article is that the prediction error would grow over many time steps
possibly leading to significant accuracy loss in the long-time dynamics. As
shown above such errors can be made small enough to make long-time predictions
reliable all the way until the dynamics reaches equilibrium.

###  IV.2 Asymmetric spin-boson model

The results reported so far are encouraging but not discriminative. To unravel
the differences between the studied ML models, we devise a more stringent
test. All the ANN and KRR models are re-trained on the asymmetric spin-boson
model data set without any hyperparameter adjustment. The increased difficulty
of this test stems from the richer dynamics of the asymmetric spin-boson model
which might require more training parameters or even different ANN
architecture (more hidden layers, longer hidden state vectors, and memory
$T$).

The proposed test clearly indicates that some models predict the long-time
dynamics much more accurately than others. The performance of each model on
the hold-out test set is reported in Tables 1 and 2. There is a noticeable
increase in the MAE of all studied ML models compared to the symmetric spin-
boson model. Generally, KRR models are more robust than ANN models. On
average, the MAE of KRR models increases by a factor of 6. KRR models based on
Mátern kernels with $n=2,3,4$ exhibit the most significant performance drop of
$\sim$10 times. The average decrease in performance of ANN models is
$\sim$18\. Interestingly, the most significant performance reduction, among
the ANN models, is observed for the models that do not contain recurrent
layers: 1D CNN and FFNN while models containing bidirectional recurrent layers
revealed to be the most robust ANN models. On average the MAE of B(LSTM, GRU,
RNN) models decreased by a factor of $\sim$10 while, for example, the FFNN
model became $\sim$30 times less accurate for the assymetric spin-boson model.

Figure 5: Expectation values $\langle\hat{\sigma}_{z}(t)\rangle$ for
asymmetric spin-boson model with $\epsilon=1$ as a function of time. Results
predicted by an indicated ML model are compared to the HEOM results (black
triangles) for the following parameters: (a)
$\lambda=0.1,\omega_{c}=6.0,\beta=0.75$; (b)
$\lambda=0.3,\omega_{c}=8.0,\beta=1.0$; (c)
$\lambda=0.2,\omega_{c}=10.0,\beta=0.25$; (d)
$\lambda=0.4,\omega_{c}=8.0,\beta=0.75$; (e)
$\lambda=0.8,\omega_{c}=10.0,\beta=1.0$; (f)
$\lambda=0.7,\omega_{c}=10.0,\beta=0.1$. All parameters are in the units of
$\Delta$.

Fig. 5 shows the ML-predicted dynamics of $\langle\hat{\sigma}_{z}(t)\rangle$
of the most and least accurate ANN and KRR models for the asymmetric spin-
boson model as well as of the “reference” ANN model, 1D CNN. Analogously to
the symmetric spin-boson model, BRNN model is found to be least accurate ANN
model. The difference is, however, more dramatic. Clearly, BRNN model is
overall no longer acceptable ANN model even though it provides an accurate
prediction for some parameters, see Fig. 5d and e, but even in these cases one
can easily distinguish an unphysical oscillatory behavoir building up beyond
relatively short times of $t\Delta\approx 10$. All other ANN models provide
more accurate and stable predictions without noise and unphysical artifacts.
The CGRU model is the most accurate matching the HEOM dynamics nearly exactly.
The worst KRR model is, once again, the one with the linear kernel. It is
however, able to capture the dynamics qualitatively, and in some, cases fairly
accurately. It seems to be less reliable in predicting the oscillatory
dynamics typically observed in low temperature and small reorganization energy
regimes. The most accurate KRR model is the one with the Gaussian kernel. This
model predicts the reference HEOM dynamics very accurately. The promise of the
KRR-G model has been already pointed out in Ref. 33.

###  IV.3 Training and prediction times

The accuracy and robustness of ML models are clearly very important but our
analysis would be incomplete without discussing the timings associated with
all studied ML models. In large-scale applications, the choice of ML model is
often a compromise between the accuracy and the complexity of the model. The
latter directly affects the training as well as the prediction times which are
important, especially for the approaches based on recursive application of ML
models to generate the dynamics. Training times for all ML models are
calculated using 2x20C Intel Xeon Gold 6230 2.1GHz processor with 192 Gb DDR4
memory, and reported in Tables 1 and 2 (only times for training are shown, in
practice, the cost is higher if the hyperparameter optimization is performed).
One should be mindful that different software packages are used for ANN and
KRR models which makes the direct comparison somewhat dubious. Nonetheless
both keras and MLatom are state-of-the-art software widely used in the
community so we will proceed with the comparison.

According to Table 1, the accuracy of FFNN and CBGRU models for the asymmetric
spin-boson model is the same but the training time of the FFNN model is nearly
4 times shorter. Furthermore, for the symmetric spin-boson Hamiltonian, FFNN
model outperforms all ANN models comprised of only recurrent layers but shows
similar performance to convolutional recurrent neural networks and
bidirectional recurrent neural networks and it does so faster: FFNN model is 3
times faster to train than the best performing CLSTM and CGRU models.

Continuing the comparison of the training times we note that, due to their
simplicity, the vanilla (unidirectional) RNN models are also fast to train.
Oppositely, bidirectional RNNs of all kinds are among the slowest models to
train. Given that using bidirectional architecture does not reduce the error
compared to their unidirectional counterparts and the increased training time
of bidirectional layers makes such models disfavoured for, at least, the
problem of predicting the long-time dynamics of the spin-boson problem. This
result is interesting given that bidirectional recurrent neural networks have
some success in other domains as described in Introduction.

Training times of the KRR models is relatively fast compared to most of the NN
models except for FFNN and CRNN.

The average single-step prediction times for each ANN and KRR models are shown
in the last column of Tables 1 and 2. It is clear that KRR models take less
time to make a single time-step prediction than ANN models, although this can
be attributed to the software nuances. Expectedly, the FFNN model is the
fastest ANN model while BGRU is the slowest. Similar to the training time we
observe that (C,B)GRU models are slower than (C,B)LSTM models which is
counerintuitive given than GRU cell is less complex than the LSTM one. We
attribute this inconsistency to the implementation details of both cells in
keras. Overall, recurrent layers require more time to evaluate than other
layers as should be anticipated. Bidirectional recurrent neural networks take
even more time to make a prediction compared to their unidirectional
counterparts. Convolutional recurrent neural network models take less time to
predict compared to their two-layer recurrent counterparts as they replace one
computationally more expensive recurrent layer with a chaper 1D CNN layer.

Figure 6: The mean absolute single-step prediction error of
$\langle\hat{\sigma}_{z}(t)\rangle$ of machine learning models studied in this
work plotted against the corresponding training (a) and single-step prediction
(b) times. Bidirectional recurrent neural network model is not shown because
its error is much too big compared to other models (see Table 1 and Fig. 5).

###  IV.4 Model efficiency

To illustrate the tradeoff between accuracy and training/prediction time, in
Fig. 6 the training and single-step prediction times are plotted against the
MAE for each ML model built in the present work. For ease of illustration, the
models of similar kind are highlighed with the same color, e.g., all KRR
models are depicted by green circles, magenta color is used for all
convolutional recurrent models, etc. Clearly, KRR models are the most
efficient models showing high accuracy and requiring the shortest amount of
time to train and make a prediction. Among the ANN models, convolutional
recurrent models are the most efficient. According to Fig. 6 convolutional
bidirectional NNs, specifically, CBLSTM constitute the second most efficient
set of ANN models. In contrast, because of the longer training times and not
much improved accuracy, bidirectional recurrent ANN models are the least
efficient.

##  V Conclusions

We performed a benchmark study of 22 supervised machine learning methods
comparing their ability to accurately forecast long-time dynamics of a two-
level quantum system linearly coupled to harmonic bath. We illustrate that
many of the studied ML methods can achieve a good agreement with the exact
population dynamics. Our study reveals that if the memory time of a quantum
dynamical system under study is known, the models based on Kernel Ridge
Regression are the most accurate and should be preferred. The commonly
employed Gaussian kernel is confirmed to be the most efficient yielding the
accuracy on a par with other nonlinear kernels while requiring somewhat less
than training time.

Convolutional recurrent neural networks appear to be the most promising ANN
models. One scenario where such models are suitable are the problems when the
memory time of the problem is not available and flexible ML models allowing
variable size input are needed. Such investigation is beyond the scope of this
work and will be performed and reported in future studies. Based on the
present study we conclude that particularly CGRU is the most promising ANN
models for long-time quantum dynamics simulations.

All ANN models are however less accurate than KRR methods with nonlinear
kernels and take more time to train and predict. Often quoted poor
computational scaling of KRR with the increasing system size may pose some
technical problems for large data sets which may be needed for large systems,
longer and larger number of training trajectories. Nevertheless, many
approaches have been suggested to mitigate this problemSnelson and Ghahramani
(2005); Deringer  _et al._ (2021); Rahimi and Recht (2007); Yu  _et al._
(2016); Hu  _et al._ (2018); Browning  _et al._ (2022) and we are also working
on implementation of alternative approaches for large data sets.

One can reasonably anticipate that our results might depend on the strategy
used to compare ML models based on artificial neural networks. Here we chosen
to build and compare ANN models with approximately the same number of
trainable parameters. There is however no proven best way to compare the
performance of ANNs with fundamentally different types of layers. We believe
that the strategy chosen in this work should provide a faithful comparison of
ANN models. Additionally, our conclusion advocating KRR methods holds
irrespective of the strategy used to compare ANN models because KRR models
studied in this work, are both faster and more accurate than ANN models.
Addition of more layers to ANN models will likely make them more accurate
(although overfitting should be carefully checked in this case) but it will
necessarily make such models more computationally expensive.

FFNN may seem as a reasonable compromise between accuracy and computational
cost. However, one needs to bear in mind that FFNN models require a fixed-size
input. Of course an input to FFNN models can be padded with zeros, but it
necessarily alters the representation of the underlying physics in the input
data and, therefore, the performance of such models is not expected to be
strong.

Many popular ML methods have been studied in the present article. However, the
field of artificial intelligence and machine learning is growing rapidly
making it nearly impossible to cover all recently developed algorithms. Our
future work will focus on novel approaches to time-series modeling such as
transformers. Vaswani  _et al._ (2017) Additionally, the CNN models employed
in this work uses kernels of the same size. Even though the kernel size was
optimized for the given task, the use of inception modules Szegedy  _et al._
(2015) that include multiple filters of varying size might improve the
performance of CNNs and will also be tested.

###### Acknowledgements.

A.A.K. acknowledges the Ralph E. Powe Junior Faculty Enhancement Award from
Oak Ridge Associated Universities. This work was also supported by the startup
funds of the College of Arts and Sciences and the Department of Physics and
Astronomy of the University of Delaware. P.O.D. acknowledges funding by the
National Natural Science Foundation of China (No. 22003051), the Fundamental
Research Funds for the Central Universities (No. 20720210092) and via the Lab
project of the State Key Laboratory of Physical Chemistry of Solid Surfaces.
K. J. R. E. acknowledges support by the Beyond Research Program between
University of Delaware and Universidad Nacional de Colombia. Calculations were
performed with high-performance computing resources provided by the University
of Delaware and Xiamen University.

## Author declarations

### Conflict of interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.

##  Appendix A Comparison of performance of KRR models trained on sub-training
and training sets

As obvious from Table 3, KRR models trained only on the sub-training set have
accuracy very close to the models trained on the entire training set (reported
in the main text). It is desirable to include all points to ensure that no
points are ”wasted”.

Table 3: Mean absolute prediction errors (MAEs), training, and average single step prediction times of all KRR models used in this work. KRR-L, KRR-G, KRR-DP, KRR-E, and KRR-M ($n=1,2,3,4$) denote kernel ridge regression models with linear kernel, Gaussian kernel, decaying-periodic kernel, exponential kernel, and Matern kernel with $n=1,2,3,4$, respectively. The MAEs are shown for both symmetric and asymmetric cases. KRR models trained on 100% training data and sub-training data (80% of training data). | Mean Absolute Prediction Error
---|---
Model | Symmetric (100 $|$ 80)% | Asymmetric (100 $|$ 80)%
KRR-L | (1.2 $|$ 1.2)$\cdot 10^{-2}$ | (6.5 $|$ 5.8)$\cdot 10^{-2}$
KRR-G | (4.7 $|$ 3.8)$\cdot 10^{-4}$ | (1.2 $|$ 1.3)$\cdot 10^{-3}$
KRR-DP | (4.3 $|$ 4.5)$\cdot 10^{-4}$ | (2.0 $|$ 2.2)$\cdot 10^{-3}$
KRR-E | (2.1 $|$ 2.5)$\cdot 10^{-3}$ | (3.3 $|$ 3.9)$\cdot 10^{-3}$
KRR-M1 | (2.4 $|$ 2.0)$\cdot 10^{-4}$ | (1.3 $|$ 1.6)$\cdot 10^{-3}$
KRR-M2 | (2.2 $|$ 3.2)$\cdot 10^{-4}$ | (2.7 $|$ 2.6)$\cdot 10^{-3}$
KRR-M3 | (2.0 $|$ 2.3)$\cdot 10^{-4}$ | (2.3 $|$ 2.1)$\cdot 10^{-3}$
KRR-M4 | (2.3 $|$ 2.4)$\cdot 10^{-4}$ | (2.1 $|$ 2.0)$\cdot 10^{-3}$

##  Appendix B Particle swarm optimization

The hyperparameter optimization of the 1D CNN model was performed using the
heuristic particle swarm optimization (PSO) algorithm. Kennedy and Eberhart
(1995) The PSO algorithm was inspired by the observation of the motion of
swarms of birds and insects, where each swarm member is guided not only by the
best solution for itself but also by the best solution seen by the entire
population.

An object of a swarm that moves around in the search space is called a
particle. A new update of the swarm is called a new generation. The particle
in the swarm includes variables (hyperparameters of the model) updated during
the optimization based on the information about the previous best (test MAE)
states of the particle and the swarm itself. Each variable of an individual
particle has velocity $\mathbf{v}(t)$, position $\mathbf{x}(t)$, and the best
particle’s position $\mathbf{x}_{p}(t)$ which has generated the smallest MAE
during the coarse of its trajectory. Additionally, there is a global variable
that contains the best global position $\mathbf{x}_{g}(t)$ of the swarm. Each
generation the position and velocity of a particle are updated according to

$\begin{array}[]{l}\mathbf{v}_{i}(t+1)=w\mathbf{v}_{i}(t)+c_{p}\left(\mathbf{x}_{p_{i}}(t)-\mathbf{x}_{i}(t)\right)+c_{g}\left(\mathbf{x}_{g_{i}}(t)-\mathbf{x}_{i}(t)\right)\\\
\mathbf{x}_{i}(t+1)=\mathbf{x}_{i}(t)+\mathbf{v}_{i}(t),\end{array}$ (39)

where $c_{p}$, $c_{g}$, and $w$ are cognitive, social, and inertia
coefficients. These coefficients quantify how much the particle is directed
toward the best solution seen by itself, by the swarm, and in the previous
direction.

In the present work, the search space is conformed by the space of
hyperparameters of the 1D CNN model as described in Sec. III.2.2. The values
of the parameters are restricted to be positive at each step. If a parameter
takes a negative value it is replaced by a random number in an interval of
values initially set for the parameter. The following values of the
coefficients are used $w=0.729$, $c_{p}=1.49445$, $c_{g}=1.49445$.

## References

  * Weiss (2012) U. Weiss,  _Quantum Dissipative Systems_ (World Scientific, Singapore, 2012).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, New York, NY, 2002).
  * Leggett  _et al._ (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, J. Phys. Soc. Jap. 58, 1199 (1989).
  * Tanimura (2020) Y. Tanimura, J. Chem. Phys. 153, 020901 (2020).
  * Meyer  _et al._ (1990) H.-D. Meyer, U. Manthe, and L. Cederbaum, Chem. Phys. lett. 165, 73 (1990).
  * Wang and Thoss (2003) H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003).
  * Makarov and Makri (1994) D. E. Makarov and N. Makri, Chem. Phys. lett. 221, 482 (1994).
  * Prior  _et al._ (2010) J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050404 (2010).
  * Ren  _et al._ (2018) J. Ren, Z. Shuai, and G. Kin-Lic Chan, J. Chem. Theory Comput. 14, 5027 (2018).
  * Luo  _et al._ (2010) B. Luo, J. Ye, C. Guan, and Y. Zhao, Phys. Chem. Chem. Phys. 12, 15073 (2010).
  * Cohen and Rabani (2011) G. Cohen and E. Rabani, Phys. Rev. B 84, 075150 (2011).
  * Cohen  _et al._ (2015) G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Phys. Rev. Lett. 115, 266802 (2015).
  * Greene and Batista (2017) S. M. Greene and V. S. Batista, J. Chem. Theory Comput. 13, 4034 (2017).
  * Yan and Shao (2016) Y.-A. Yan and J. Shao, Front. Phys. 11, 1 (2016).
  * Hsieh and Cao (2018a) C.-Y. Hsieh and J. Cao, J. Chem. Phys. 148, 014103 (2018a).
  * Hsieh and Cao (2018b) C.-Y. Hsieh and J. Cao, J. Chem. Phys. 148, 014104 (2018b).
  * Han  _et al._ (2020) L. Han, A. Ullah, Y.-A. Yan, X. Zheng, Y. Yan, and V. Chernyak, J. Chem. Phys. 152, 204105 (2020).
  * Ullah  _et al._ (2020) A. Ullah, L. Han, Y.-A. Yan, X. Zheng, Y. Yan, and V. Chernyak, J. Chem. Phys. 152, 204106 (2020).
  * Nakajima (1958) S. Nakajima, Prog. Theor. Phys. 20, 948 (1958).
  * Zwanzig (1960) R. Zwanzig, J. Chem. Phys. 33, 1338 (1960).
  * Shi and Geva (2003) Q. Shi and E. Geva, J. Chem. Phys. 119, 12063 (2003).
  * Kelly and Markland (2013) A. Kelly and T. E. Markland, J. Chem. Phys. 139, 014104 (2013).
  * Mulvihill and Geva (2021) E. Mulvihill and E. Geva, J. Phys. Chem. B 125, 9834 (2021).
  * Mulvihill  _et al._ (2021) E. Mulvihill, K. M. Lenn, X. Gao, A. Schubert, B. D. Dunietz, and E. Geva, J. Chem. Phys. 154, 204109 (2021).
  * Brian and Sun (2021) D. Brian and X. Sun, Chinese J. Chem. Phys. 34, 497 (2021).
  * Cerrillo and Cao (2014) J. Cerrillo and J. Cao, Phys. Rev. Lett. 112, 110401 (2014).
  * Kananenka  _et al._ (2016) A. A. Kananenka, C.-Y. Hsieh, J. Cao, and E. Geva, J. Phys. Chem. Lett. 7, 4809 (2016).
  * Buser  _et al._ (2017) M. Buser, J. Cerrillo, G. Schaller, and J. Cao, Phys. Rev. A 96, 062122 (2017).
  * Gelzinis  _et al._ (2017) A. Gelzinis, E. Rybakovas, and L. Valkunas, J. Chem. Phys. 147, 234108 (2017).
  * Chen  _et al._ (2020) Y.-Q. Chen, K.-L. Ma, Y.-C. Zheng, J. Allcock, S. Zhang, and C.-Y. Hsieh, Phys. Rev. Applied 13, 034045 (2020).
  * Herrera Rodríguez and Kananenka (2021) L. E. Herrera Rodríguez and A. A. Kananenka, J. Phys. Chem. Lett. 12, 2476 (2021).
  * Ullah and Dral (2021) A. Ullah and P. O. Dral, New J. Phys. 23, 113019 (2021).
  * Ullah and Dral (2022a) A. Ullah and P. O. Dral, Nature Commun. 13, 1 (2022a).
  * Ullah and Dral (2022b) A. Ullah and P. O. Dral, J. Phys. Chem. Lett. 13, 6037 (2022b).
  * Box  _et al._ (2015) G. Box, G. Jenkins, G. Reinsel, and G. Ljung,  _Time Series Analysis: Forecasting and Control_ , Wiley Series in Probability and Statistics (Wiley, 2015).
  * Ariyo  _et al._ (2014) A. A. Ariyo, A. O. Adewumi, and C. K. Ayo, in  _2014 UKSim-AMSS 16th International Conference on Computer Modelling and Simulation_ (2014) pp. 106–112.
  * Khashei and Bijari (2011) M. Khashei and M. Bijari, Appl. Soft Comput. 11, 2664 (2011).
  * Schmidhuber (2015) J. Schmidhuber, Neural Netw. 61, 85 (2015).
  * Goodfellow  _et al._ (2016) I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_, Adaptive Computation and Machine Learning series (MIT Press, 2016).
  * Sherstinsky (2020) A. Sherstinsky, Phys. D: Nonlinear Phenom. 404, 132306 (2020).
  * Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber, Neural Comput. 9, 1735 (1997).
  * Ye  _et al._ (2012) Q. Ye, W. Y. Szeto, and S. C. Wong, IEEE Transactions on Intelligent Transportation Systems 13, 1727 (2012).
  * Siami-Namini  _et al._ (2018) S. Siami-Namini, N. Tavakoli, and A. Siami Namin, in  _2018 17th IEEE International Conference on Machine Learning and Applications (ICMLA)_ (2018) pp. 1394–1401.
  * Siami-Namini  _et al._ (2019) S. Siami-Namini, N. Tavakoli, and A. S. Namin, in  _2019 IEEE International Conference on Big Data (Big Data)_ (2019) pp. 3285–3292.
  * Malhotra  _et al._ (2017) P. Malhotra, V. TV, L. Vig, P. Agarwal, and G. Shroff, arXiv preprint arXiv:1706.08838 (2017).
  * Fawaz  _et al._ (2019) H. I. Fawaz, G. Forestier, J. Weber, L. Idoumghar, and P.-A. Muller, Data Min. Knowl. Discov. 33, 917 (2019), 1809.04356 .
  * Kashiparekh  _et al._ (2019) K. Kashiparekh, J. Narwariya, P. Malhotra, L. Vig, and G. Shroff, in  _2019 International Joint Conference on Neural Networks (IJCNN)_ (IEEE, 2019) pp. 1–8.
  * Moghar and Hamiche (2020) A. Moghar and M. Hamiche, Procedia Comput. Sci. 170, 1168 (2020), the 11th International Conference on Ambient Systems, Networks and Technologies (ANT) / The 3rd International Conference on Emerging Data and Industry 4.0 (EDI40) / Affiliated Workshops.
  * Fischer and Krauss (2018) T. Fischer and C. Krauss, Eur. J. Oper. Res. 270, 654 (2018).
  * Kim and Moon (2019) J. Kim and N. Moon, J. Ambient Intell. Humaniz. Comput. , 1 (2019).
  * Kumar  _et al._ (2021) D. Kumar, P. K. Sarangi, and R. Verma, Mater. Today: Proc. (2021).
  * Graves (2012a) A. Graves,  _Supervised Sequence Labelling with Recurrent Neural Networks_, Studies in Computational Intelligence (Springer Berlin Heidelberg, 2012).
  * Graves and Schmidhuber (2005) A. Graves and J. Schmidhuber, Neural Netw. 18, 602 (2005), iJCNN 2005.
  * Graves (2012b) A. Graves, arXiv preprint arXiv:1211.3711 (2012b).
  * Graves  _et al._ (2013) A. Graves, A.-r. Mohamed, and G. Hinton, in  _2013 IEEE international conference on acoustics, speech and signal processing_ (Ieee, 2013) pp. 6645–6649.
  * Hannun  _et al._ (2014) A. Hannun, C. Case, J. Casper, B. Catanzaro, G. Diamos, E. Elsen, R. Prenger, S. Satheesh, S. Sengupta, A. Coates, _et al._ , arXiv preprint arXiv:1412.5567 (2014).
  * Xiong  _et al._ (2016) W. Xiong, J. Droppo, X. Huang, F. Seide, M. Seltzer, A. Stolcke, D. Yu, and G. Zweig, arXiv preprint arXiv:1610.05256 (2016).
  * Graves  _et al._ (2007) A. Graves, S. Fernández, M. Liwicki, H. Bunke, and J. Schmidhuber, in  _Proceedings of the 20th International Conference on Neural Information Processing Systems_ , NIPS’07 (Curran Associates Inc., Red Hook, NY, USA, 2007) p. 577–584.
  * Graves and Schmidhuber (2009) A. Graves and J. Schmidhuber, in  _Advances in Neural Information Processing Systems_, Vol. 21, edited by D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou (Curran Associates, Inc., 2009).
  * Tang  _et al._ (2016) Y. Tang, Y. Huang, Z. Wu, H. Meng, M. Xu, and L. Cai, in  _2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ (2016) pp. 6125–6129.
  * Yin  _et al._ (2017) W. Yin, K. Kann, M. Yu, and H. Schütze, arXiv preprint arXiv:1702.01923 (2017).
  * Sutskever  _et al._ (2014) I. Sutskever, O. Vinyals, and Q. V. Le, Adv. Neural Inf. Process. Syst. 27 (2014).
  * Bahdanau  _et al._ (2014) D. Bahdanau, K. Cho, and Y. Bengio, arXiv preprint arXiv:1409.0473 (2014).
  * Shido  _et al._ (2019) Y. Shido, Y. Kobayashi, A. Yamamoto, A. Miyamoto, and T. Matsumura, in  _2019 International Joint Conference on Neural Networks (IJCNN)_ (2019) pp. 1–8.
  * Choi  _et al._ (2017) E. Choi, M. T. Bahadori, L. Song, W. F. Stewart, and J. Sun, arXiv preprint arXiv:1611.07012v3 (2017).
  * Lynn  _et al._ (2019) H. M. Lynn, S. B. Pan, and P. Kim, IEEE Access 7, 145395 (2019).
  * Gupta  _et al._ (2021) P. Gupta, P. Malhotra, J. Narwariya, L. Vig, and G. Shroff, arXiv preprint arXiv:1904.00655 (2021).
  * Ma  _et al._ (2017) F. Ma, R. Chitta, J. Zhou, Q. You, T. Sun, and J. Gao, in  _Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining_ (2017) pp. 1903–1911.
  * van Lint  _et al._ (2005) J. van Lint, S. Hoogendoorn, and H. van Zuylen, Transp. Res. Part C Emerg. Technol. 13, 347 (2005).
  * Zhao  _et al._ (2017) Z. Zhao, W. Chen, X. Wu, P. C. Y. Chen, and J. Liu, IET Intell. Transp. Syst. 11, 68 (2017).
  * Cui  _et al._ (2019) Z. Cui, R. Ke, Z. Pu, and Y. Wang, arXiv preprint arXiv:1801.02143 (2019).
  * Eck and Schmidhuber (2002) D. Eck and J. Schmidhuber, “A first look at music composition using lstm recurrent neural networks,” (2002), istituto Dalle Molle Di Studi Sull Intelligenza Artificiale.
  * Boulanger-Lewandowski  _et al._ (2012) N. Boulanger-Lewandowski, Y. Bengio, and P. Vincent, arXiv preprint arXiv:1206.6392 (2012).
  * Srivastava  _et al._ (2015) N. Srivastava, E. Mansimov, and R. Salakhudinov, in  _Proceedings of the 32nd International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 37, edited by F. Bach and D. Blei (PMLR, Lille, France, 2015) pp. 843–852.
  * Habi and Messer (2019) H. V. Habi and H. Messer, in  _2019 IEEE International Workshop on Signal Processing Systems (SiPS)_ (2019) pp. 184–188.
  * Gupta  _et al._ (2018) A. Gupta, A. T. Müller, B. J. H. Huisman, J. A. Fuchs, P. Schneider, and G. Schneider, Mol. Inform. 37, 1700111 (2018).
  * Segler  _et al._ (2018) M. H. S. Segler, T. Kogej, C. Tyrchan, and M. P. Waller, ACS Cent. Sci. 4, 120 (2018).
  * Wei and Yan (2019) W. Wei and X. Yan, IOP Conf. Ser.: Mater. Sci. Eng. 688, 033022 (2019).
  * Abbasimehr  _et al._ (2020) H. Abbasimehr, M. Shabani, and M. Yousefi, Comput. Ind. Eng. 143, 106435 (2020).
  * Zhang  _et al._ (2020) L. Zhang, C. Xu, Y. Gao, Y. Han, X. Du, and Z. Tian, Tsinghua Sci. Technol. 25, 712 (2020).
  * Lyu  _et al._ (2016) H. Lyu, H. Lu, and L. Mou, Remote Sens. 8, 506 (2016).
  * Ienco  _et al._ (2017) D. Ienco, R. Gaetano, C. Dupaquier, and P. Maurel, IEEE Geoscience and Remote Sensing Letters 14, 1685 (2017).
  * Bhoopchand  _et al._ (2016) A. Bhoopchand, T. Rocktäschel, E. Barr, and S. Riedel, arXiv preprint arXiv:1611.08307 (2016).
  * Pineda (1987) F. J. Pineda, Phys. Rev. Lett. 59, 2229 (1987).
  * Gers  _et al._ (2000) F. A. Gers, J. Schmidhuber, and F. Cummins, Neural Comput. 12, 2451 (2000).
  * Gers  _et al._ (1999) F. Gers, J. Schmidhuber, and F. Cummins, in  _1999 Ninth International Conference on Artificial Neural Networks ICANN 99. (Conf. Publ. No. 470)_, Vol. 2 (1999) pp. 850–855 vol.2.
  * Cho  _et al._ (2014) K. Cho, B. Van Merriënboer, D. Bahdanau, and Y. Bengio, arXiv preprint arXiv:1409.1259 (2014).
  * Chung  _et al._ (2016) J. Chung, K. Cho, and Y. Bengio, arXiv preprint arXiv:1603.06147 (2016).
  * Jozefowicz  _et al._ (2015) R. Jozefowicz, W. Zaremba, and I. Sutskever, in  _Proceedings of the 32nd International Conference on Machine Learning_ , Proceedings of Machine Learning Research, Vol. 37, edited by F. Bach and D. Blei (PMLR, Lille, France, 2015) pp. 2342–2350.
  * Schuster and Paliwal (1997) M. Schuster and K. Paliwal, IEEE Trans. Signal Process. 45, 2673 (1997).
  * Baldi  _et al._ (1999) P. Baldi, S. Brunak, P. Frasconi, G. Soda, and G. Pollastri, Bioinformatics 15, 937 (1999).
  * LeCun  _et al._ (1989) Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel, Neural Comput. 1, 541 (1989).
  * Wang  _et al._ (2016) Z. Wang, W. Yan, and T. Oates, arXiv preprint arXiv:1611.06455 (2016).
  * Cui  _et al._ (2016) Z. Cui, W. Chen, and Y. Chen, arXiv preprint arXiv:1603.06995 (2016).
  * Serrà  _et al._ (2018) J. Serrà, S. Pascual, and A. Karatzoglou, arXiv preprint arXiv:1805.03908 (2018).
  * Zheng  _et al._ (2014) Y. Zheng, Q. Liu, E. Chen, Y. Ge, and J. L. Zhao, in  _Web-Age Information Management_ , edited by F. Li, G. Li, S.-w. Hwang, B. Yao, and Z. Zhang (Springer International Publishing, Cham, 2014) pp. 298–310.
  * Zheng  _et al._ (2016) Y. Zheng, Q. Liu, E. Chen, Y. Ge, and J. L. Zhao, Front. Comput. Sci. 10, 96 (2016).
  * Interdonato  _et al._ (2018) R. Interdonato, D. Ienco, R. Gaetano, and K. Ose, arXiv preprint arXiv:1809.07589 (2018).
  * Ismail Fawaz  _et al._ (2020) H. Ismail Fawaz, B. Lucas, G. Forestier, C. Pelletier, D. F. Schmidt, J. Weber, G. I. Webb, L. Idoumghar, P.-A. Muller, and F. Petitjean, Data Min Knowl Discov 34, 1936–1962 (2020).
  * Tang  _et al._ (2021) W. Tang, G. Long, L. Liu, T. Zhou, J. Jiang, and M. Blumenstein, arXiv preprint arXiv:2002.10061 (2021).
  * Rajpurkar  _et al._ (2017) P. Rajpurkar, A. Y. Hannun, M. Haghpanahi, C. Bourn, and A. Y. Ng, arXiv preprint arXiv:1707.01836 (2017).
  * Roy  _et al._ (2018) S. Roy, I. Kiral-Kornek, and S. Harrer, arXiv preprint arXiv:1802.00308 (2018).
  * Schirrmeister  _et al._ (2018) R. T. Schirrmeister, L. Gemein, K. Eggensperger, F. Hutter, and T. Ball, arXiv preprint arXiv:1708.08012 (2018).
  * Sercu  _et al._ (2016) T. Sercu, C. Puhrsch, B. Kingsbury, and Y. LeCun, arXiv preprint arXiv:1509.08967 (2016).
  * Choi  _et al._ (2016) K. Choi, G. Fazekas, M. Sandler, and K. Cho, arXiv preprint arXiv:1609.04243 (2016).
  * You  _et al._ (2018) S. D. You, C.-H. Liu, and W.-K. Chen, Hum.-centric Comput. Inf. Sci. 8, 34 (2018).
  * Karim  _et al._ (2018) F. Karim, S. Majumdar, H. Darabi, and S. Chen, IEEE Access 6, 1662–1669 (2018).
  * Sheykhivand  _et al._ (2020) S. Sheykhivand, Z. Mousavi, T. Y. Rezaii, and A. Farzamnia, IEEE Access 8, 139332 (2020).
  * Eapen  _et al._ (2019) J. Eapen, D. Bein, and A. Verma, in  _2019 IEEE 9th Annual Computing and Communication Workshop and Conference (CCWC)_ (2019) pp. 0264–0270.
  * Shi  _et al._ (2015) X. Shi, Z. Chen, H. Wang, D.-Y. Yeung, W. kin Wong, and W. chun Woo, arXiv preprint arXiv:1506.04214 (2015).
  * Müller  _et al._ (1997) K. R. Müller, A. J. Smola, G. Rätsch, B. Schölkopf, J. Kohlmorgen, and V. Vapnik, in  _Artificial Neural Networks — ICANN’97_ , edited by W. Gerstner, A. Germond, M. Hasler, and J.-D. Nicoud (Springer Berlin Heidelberg, Berlin, Heidelberg, 1997) pp. 999–1004.
  * Sapankevych and Sankar (2009) N. I. Sapankevych and R. Sankar, IEEE Comput. Intell. Mag. 4, 24 (2009).
  * Haworth  _et al._ (2014) J. Haworth, J. Shawe-Taylor, T. Cheng, and J. Wang, Transp. Res. C: Emerg. Technol. 46, 151 (2014).
  * Gneiting  _et al._ (2010) T. Gneiting, W. Kleiber, and M. Schlather, J. Am. Stat. Assoc. 105, 1167 (2010).
  * Cuturi and Doucet (2011) M. Cuturi and A. Doucet, arXiv preprint arXiv:1101.0673 (2011).
  * Akimov (2021) A. V. Akimov, J. Phys. Chem. Lett. 12, 12119 (2021).
  * Secor  _et al._ (2021) M. Secor, A. V. Soudackov, and S. Hammes-Schiffer, J. Phys. Chem. Lett. 12, 10654 (2021).
  * Yang  _et al._ (2020) B. Yang, B. He, J. Wan, S. Kubal, and Y. Zhao, Chem. Phys. 528, 110509 (2020).
  * Bandyopadhyay  _et al._ (2018) S. Bandyopadhyay, Z. Huang, K. Sun, and Y. Zhao, Chem. Phys. 515, 272 (2018).
  * Banchi  _et al._ (2018) L. Banchi, E. Grant, A. Rocchetto, and S. Severini, New J. Phys. 20, 123030 (2018).
  * Wu  _et al._ (2021) D. Wu, Z. Hu, J. Li, and X. Sun, J. Chem. Phys. 155, 224104 (2021).
  * Lin  _et al._ (2021) K. Lin, J. Peng, F. L. Gu, and Z. Lan, J. Phys. Chem. Lett. 12, 10225 (2021).
  * Choi  _et al._ (2022) M. Choi, D. Flam-Shepherd, T. H. Kyaw, and A. Aspuru-Guzik, Phys. Rev. A 105, 042403 (2022).
  * Tsai  _et al._ (2022) S.-T. Tsai, E. Fields, Y. Xu, E.-J. Kuo, and P. Tiwary, arXiv preprint arXiv:2203.00597 (2022).
  * Tang  _et al._ (2022) D. Tang, L. Jia, L. Shen, and W.-H. Fang, arXiv preprint arXiv:2206.13780 (2022).
  * Dral and Barbatti (2021) P. O. Dral and M. Barbatti, Nat. Rev. Chem. 5, 388–405 (2021).
  * Westermayr and Marquetand (2020) J. Westermayr and P. Marquetand, Mach. Learn.: Sci. Technol. 1, 043001 (2020).
  * Westermayr and Marquetand (2021) J. Westermayr and P. Marquetand, Chem. Rev. 121, 9873–9926 (2021).
  * Lin  _et al._ (2022) K. Lin, J. Peng, C. Xu, F. L. Gu, and Z. Lan, arXiv preprint arXiv:2205.03600 (2022).
  * Makhlin  _et al._ (2001) Y. Makhlin, G. Schön, and A. Shnirman, Rev. Mod. Phys. 73, 357 (2001).
  * Alvermann and Fehske (2009) A. Alvermann and H. Fehske, Phys. Rev. Lett. 102, 150601 (2009).
  * Winter  _et al._ (2009) A. Winter, H. Rieger, M. Vojta, and R. Bulla, Phys. Rev. Lett. 102, 030601 (2009).
  * Garg  _et al._ (1985) A. Garg, J. N. Onuchic, and V. Ambegaokar, J. Chem. Phys. 83, 4491 (1985).
  * Wang  _et al._ (1999) H. Wang, X. Song, D. Chandler, and W. H. Miller, J. Chem. Phys. 110, 4828 (1999).
  * Hornik  _et al._ (1989) K. Hornik, M. Stinchcombe, and H. White, Neural Netw. 2, 359 (1989).
  * Cybenko (1989) G. Cybenko, Math. Control Signals Syst. 2, 303 (1989).
  * Leshno  _et al._ (1993) M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, Neural Netw. 6, 861 (1993).
  * Zhou and Chellappa (1988) Zhou and Chellappa, in  _IEEE 1988 International Conference on Neural Networks_ (1988) pp. 71–78 vol.2.
  * Rumelhart  _et al._ (1986) D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Nature 323, 533 (1986).
  * Schäfer and Zimmermann (2006) A. M. Schäfer and H. G. Zimmermann, in  _Artificial Neural Networks – ICANN 2006_ , edited by S. D. Kollias, A. Stafylopatis, W. Duch, and E. Oja (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 632–640.
  * Chollet  _et al._ (2015) F. Chollet _et al._ , “Keras,” (2015), https://github.com/fchollet/keras.
  * Hochreiter (1991) J. Hochreiter, “Untersuchungen zu dynamischen neuronalen netzen,” Ph.D. thesis (1991), technische Universitat Munchen.
  * Bengio  _et al._ (1994) Y. Bengio, P. Simard, and P. Frasconi, IEEE trans. neural netw. 5, 157 (1994).
  * Bengio  _et al._ (1993) Y. Bengio, P. Frasconi, and P. Y. Simard, in  _Proceedings of International Conference on Neural Networks (ICNN’88), San Francisco, CA, USA, March 28 - April 1, 1993_ (IEEE, 1993) pp. 1183–1188.
  * Pascanu  _et al._ (2013) R. Pascanu, T. Mikolov, and Y. Bengio, arXiv preprint arXiv:1211.5063 (2013).
  * Sussillo and Abbott (2014) D. Sussillo and L. Abbott, arXiv preprint arXiv:1412.6558 (2014).
  * Chung  _et al._ (2015) J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, in  _International conference on machine learning_ (PMLR, 2015) pp. 2067–2075.
  * Abadi  _et al._ (2015) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,”  (2015), software available from tensorflow.org.
  * Stulp and Sigaud (2015) F. Stulp and O. Sigaud, Neural Netw. 69, 60 (2015).
  * Dral (2019a) P. O. Dral, J. Comput. Chem. 40, 2339 (2019a).
  * Hastie  _et al._ (2009) T. Hastie, R. Tibshirani, J. H. Friedman, and J. H. Friedman,  _The elements of statistical learning: data mining, inference, and prediction_ , Vol. 2 (Springer, 2009).
  * Rasmussen and Williams (2005) C. Rasmussen and C. Williams,  _Gaussian Processes for Machine Learning_, Adaptive Computation and Machine Learning series (MIT Press, 2005).
  * Dral (2019b) P. O. Dral, J. Comput. Chem. 40, 2339 (2019b).
  * Dral  _et al._ (2021) P. O. Dral, F. Ge, B.-X. Xue, Y.-F. Hou, M. Pinheiro, J. Huang, and M. Barbatti, Top. Curr. Chem. 379, 1 (2021).
  * Dral  _et al._ (2022) P. O. Dral, P. Zheng, B.-X. Xue, F. Ge, Y.-F. Hou, and M. Pinheiro Jr,  _MLatom: a Package for Atomistic Simulations with Machine Learning, development version.http://MLatom.com (accessed on June 29, 2022)_, Xiamen University, Xiamen 361005, China (2013–2022).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, _et al._ , J. Mach. Learn. Res. 12, 2825 (2011).
  * Johansson  _et al._ (2012) J. Johansson, P. Nation, and F. Nori, Comput. Phys. Commun. 183, 1760 (2012).
  * Ullah and Dral (2022c) A. Ullah and P. O. Dral, Nature Commun. 13, 1930 (2022c).
  * Kingma and Ba (2017) D. P. Kingma and J. Ba, arXiv preprint arXiv:1412.6980 (2017).
  * Glorot and Bengio (2010) X. Glorot and Y. Bengio, in  _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_, Proceedings of Machine Learning Research, Vol. 9, edited by Y. W. Teh and M. Titterington (PMLR, Chia Laguna Resort, Sardinia, Italy, 2010) pp. 249–256.
  * Bergstra  _et al._ (2011) J. Bergstra, R. Bardenet, Y. Bengio, and B. Kégl, “Algorithms for hyper-parameter optimization,” in  _Advances in Neural Information Processing Systems_, Vol. 24, edited by J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Q. Weinberger (Curran Associates, Inc., 2011).
  * Bergstra  _et al._ (2015) J. Bergstra, B. Komer, C. Eliasmith, D. L. K. Yamins, and D. D. Cox, Comput. Sci. Discov. 8, 014008 (2015).