Machine learning classification of non-Markovian noise disturbing quantum dynamics

For the generation of the data used to train the ML-models, we consider two variants of three different classification problems. Each sample of the data sets is created by first generating a random set of links ${ \mathcal E }$ (random topology) for the graph ${ \mathcal G }$, and then initialising the particle in a randomly chosen node of the graph. We set M = 15 as the number of evaluations (measurements) of the quantum particle dynamics, and d = 40 as the number of nodes of the graph ${ \mathcal G }$. This means that ${{ \mathcal P }}_{{t}_{0}}$ is a Kronecker delta centered in one of the 40 nodes, and the stochastic quantum dynamics is evolved for 15 steps for each simulated noise source of the generated data set. Here, it is worth noting that the choice of M = 15 is dictated by the fact that in recent experiments as for instance in [9, 54, 55], the number of intermediate quantum measurements does not exceed 10, and thus M = 15 is sufficiently large to represent actual physical setups. Instead, regarding taking d = 40, such a value is just able to generate a complex landscape for the particle dynamics and small enough to be numerically manageable. The total considered dynamical time t_M is taken equal to t_M = 1 or t_M = 0.1 in dimensionless units, each of them corresponding to a specific variant. Notice, indeed, that the values of t_M are expressed consistently with the energy scale of the couplings g_t , whose random values g^(j) belong to the set {1, 2, 3, 4, 5} in the data set generation, such that ℏ can be reliably set to 1 as usual. All the probability distributions ${{ \mathcal P }}_{{t}_{k}}$ for k = 0, ... , 15 are stored together with the attached label that indicates the associated type of noise.

For each of the two variants of our classification problems, we generate three different balanced data sets of 20 000 samples. The first data set, which we call IID, is suitable for a supervised binary classification task that discriminates between two different i.i.d. noisy quantum dynamics, where the noise sources have the same support but different probability distribution Prob(g).

The second data set, named as NM, concerns the classification of two different coloured noisy quantum dynamics with the noise sources again having the same support but different Prob(g) (the same ones as in the data set IID) and a transition matrix T.

Finally, the third data set, called VS, is created for the classification between stochastic quantum dynamics affected respectively by an i.i.d. and a coloured noise with same support and Prob(g).

Note that choosing graphs with random links allows to increase the statistical variability of the input data, with the result that the ML algorithms learn to classify noise sources independently of the graph topology. The aim, indeed, is to prevent that the ML-models rely only on features specific to a small class of topologies. Moreover, taking random initial distributions ${{ \mathcal P }}_{{t}_{0}}$ allows to increase the robustness of the ML methods, making them less likely to overfit on the synthetic data set.

As it will be explained in the following, some ML-models that we are going to introduce will use as input only the last distribution ${{ \mathcal P }}_{{t}_{15}}$, while other ML-models will take all the ${{ \mathcal P }}_{{t}_{k}}$ for any t_k . Moreover, each data set is balanced split in a training set of 12 000 samples, a validation set of 4 000 samples, and a test set of 4 000 samples.

In tables 1 and 2 we plot the occupation probabilities ${{ \mathcal P }}_{{t}_{k}}$ (just for the IID case for the sake of an easier presentation), being here interested in looking for the difference between choosing t₁₅ = 0.1 or 1, which identify the two different variants of the generated data set. In this regard, it is worth noting that the duration t₁₅ = 0.1 (in dimensionless units) of the quantum system dynamics, as in the example in table 1, is the minimal one to observe the diffusion of the system's population outside the node on which has been initialised. However, as it will be verified by our experiments and explained later, with this choice one has that, by taking t₁₅ = 0.1, the classification problem results quite straightforward. Indeed, just basic ML-models that are only trained on ${{ \mathcal P }}_{{t}_{15}}$ (thus, only on the final distribution ${ \mathcal P }$) are able to correctly classify between two noisy quantum dynamics. Therefore, it was more interesting to increase the value of t₁₅ up to t₁₅ = 1 (in dimensionless units). As in the example of table 2, it leads to more complex data sets, and only deep learning models, designed to read all the ${{ \mathcal P }}_{{t}_{k}}$, can classify the generated noisy quantum dynamics.

Table 1. Example of a part of ${{ \mathcal P }}_{{t}_{k}}$ for all the discrete time instants t_k for a noisy quantum dynamics affected by i.i.d. noise sources and t₁₅ = 0.1 (in dimensionless units). In the table, ${{ \mathcal P }}_{{t}_{k}}^{(s)}$ denotes the s-th element of the vector ${{ \mathcal P }}_{{t}_{k}}$ for any t_k , k = 0, ... , 15.

Table 2. Example of a part of ${{ \mathcal P }}_{{t}_{k}}$ for all the discrete time instants t_k for a noisy quantum dynamics affected by i.i.d. noise sources and t₁₅ = 1 (in dimensionless units). Again, ${{ \mathcal P }}_{{t}_{k}}^{(s)}$ denotes the s-th element of the vector ${{ \mathcal P }}_{{t}_{k}}$ for any t_k , k = 0, ... , 15. The topology and the initial state, for this example, are the same of those in table 1.

As final remark, note that the current synthetic data set is build assuming perfect measurement statistics, as it was obtained from a large enough number of repetitions of the noisy quantum dynamics. Hence, to better adapt the synthetic data set to real data, one should simulate experimental case in which the measurement statistics are estimated from a finite number of dynamics realizations (i.e., measurement shots).

We here present the binary supervised classification tasks that we are going to address, by taking ${{ \mathcal P }}_{{t}_{k}}$ as input:

• Two different probability distributions Prob(g)—specifically, ${p}_{g}^{(1)},\,\ldots \,,{p}_{g}^{(5)}=(0.0124,0.04236,0.0820\,,0.2398,0.6234)$ and =(0.1782, 0.1865, 0.2, 0.2107, 0.2245)—both associated with i.i.d. noise sources.
• Two different Prob(g) (the same as (i)) and different values of the correlation parameters—identified by transition matrices T as explained in section 2—for coloured (thus, non-Markovian) noise processes.
• An i.i.d. and a coloured noise process with the same support g and distribution ${p}_{g}^{(1)},\,\ldots \,,{p}_{g}^{(5)}=(0.0124\,,0.04236,0.0820,0.2398,0.6234)$ that thus differ for the presence of non-zero correlation parameters.

The values of both Prob(g) and the transition matrices T, used in our numerical simulations, are chosen randomly.

To solve the classification tasks above, we have employed in this paper both a more standardized ML-model that is Support Vector Machine (SVM) [33] and more recent Artificial Neural Networks (ANNs) [47, 56–58] models in the form of MLPs and RNNs. For an exhaustive explanation of such a ML-models refer to the appendix.

Let us now investigate the scaling of the classification accuracy γ, as a function of both the interval Δ between two consecutive transitions for g and the number M of discrete time instants. Notice that Δ and M are related to the total dynamical time t_M , since t_M ≡ MΔ.

A possible explanation of the differences observed between the three previously-analysed scenarios, i.e., t₁₅ = 0.1, t₁₅ = 1 and t₁₅ ≫ 1 (in dimensionless units), could be that the information on both the noise source and the initial quantum state is lost during the evolution of the system. For such aspect, not only the total dynamical time t₁₅ could play a role, but also the time interval Δ ≡ t₁ − t₀ ≡ ⋯ ≡ t_M − t_M−1. In fact, it is reasonable to conjecture that a ML-model, able to correctly classify our noisy quantum dynamics with t₁₅ = 1 (thus M = 15), can also work with ${t}_{M^{\prime} }\gg 1$ for $M^{\prime} \gt 15$ and ${\rm{\Delta }}^{\prime} ={\rm{\Delta }}$ where ${\rm{\Delta }}^{\prime} \equiv {t}_{1}-{t}_{0}\equiv \cdots \equiv \,{t}_{M^{\prime} }-{t}_{M^{\prime} -1}$. In this way, the sequence ${{ \mathcal P }}_{{t}_{1}},\,\ldots \,,{{ \mathcal P }}_{{t}_{15}}$ is contained in ${{ \mathcal P }}_{{t}_{1}},\,\ldots \,,{{ \mathcal P }}_{{t}_{M^{\prime} }}$. In other terms, we conjecture that the classification problem can be solved even for longer noisy quantum dynamics, but provided that Δ remains small.

To gain evidence on this conjecture, we have performed two additional experiments. Starting from the task IID with t₁₅ = 1 and m-biGRU-max as baseline (accuracy 91.8%), the same model (optimised in the same hyperparameters space) is trained on two new data sets. In both data sets, t_M = 2 with M equal to 15 for the first data set and 30 for the second one. Thus, in the former ${\rm{\Delta }}^{\prime} \gt {\rm{\Delta }}$ with Δ time interval of the original data set, while in the latter ${\rm{\Delta }}^{\prime} ={\rm{\Delta }}$. The first experiment (${\rm{\Delta }}^{\prime} \gt {\rm{\Delta }}$) provides a classification accuracy of 81.1%, contrarily to the results from the second experiment (${\rm{\Delta }}^{\prime} ={\rm{\Delta }}$), where a better accuracy of 96.3% is achieved. We thus observe that, by taking ${\rm{\Delta }}^{\prime} ={\rm{\Delta }}$ and the same ML-model, the classification problem can be solved with an higher accuracy, but at the price of a longer training time. Indeed, in this case, the length of each sample of the data set is twice the original one.

In another experiment, whose results are shown in figure 2, we vary M by keeping the total evolution time equal to t_M = 1. Such tests use as baseline the model m-biGRU-max applied on the most difficult task of table 3, i.e., NM with t₁₅ = 1. As a result, the achieved classification accuracy is directly proportional to M and, thus, inversely proportional to the value of the time interval Δ. Indeed, by taking t_M fixed and reducing Δ, the classification accuracy of the same model can be enhanced. Specifically, it is possible to obtain more than 90% of accuracy also for the task NM with a total dynamical time equal to 1, at the price of a longer training time as the length of the sequences increases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here, we address the following question: Could the proposed ML techniques be applied for the inference of noise sources affecting the dynamics of classical systems, e.g., Langevin equations [59]? Probably yes, but we expect that their application to quantum systems, maybe surprisingly, can be more effective than on classical systems. Both classical (non-periodic) dissipative dynamics [59] and stochastic quantum dynamics (stochastic due to the presence of an external environment, or noise sources as in our case) can asymptotically tend to a fixed-point, whereby the information on the initial state is lost. This means that the states of the system used for this noise classification tend to become indistinguishable as time increases. Classically, this can happen due to energy dissipation introduced by damping terms Instead, quantum-mechanically, a dynamical fixed point can be reached due to decoherence that makes vanishing, at least on average, all quantum coherence terms [5]. Thus, once the transient of the evolution is elapsed, the evaluation of the final state of the system does not bring information neither on the initial state nor on the initial dynamics bringing the system to the asymptotic fixed-point. In our case, we have observed that, by using only ${{ \mathcal P }}_{{t}_{15}}$ with long total dynamical time, the accuracy of all the classification tasks is always around 50% both for SVM and MLP. Consequently, if one aims to infer/reconstruct the value of parameters, signals or operators that influence the system dynamics by measuring its evolution, the most appropriate time window is during the transient. In this regard, a quantum dynamic, until it is nearly close of being unitary, is able to explore different configurations thanks to linearity and the quantum superposition principle. Conversely, classical dynamics, not being able to propagate superpositions of their trajectories, cannot provide per time unit the same amount of information on the quantity to be inferred.

In conclusion, the application of the proposed methods should be more accurate if applied to quantum systems than classical ones, but during the transient of its dynamical evolution when quantum effects are still predominant and the distance among the state and the fixed point is not negligible.

Our techniques are expected to be adopted to witness the (non-)Markovianity of the noise sources in commercial quantum devices, as for example the Q-IBM® [49] or Rigetti®. In fact, such quantum devices, as other Noisy Intermediate-Scale Quantum prototypes [60], are unavoidably affected by the external environment that entails random errors. Recently, in [61, 62], it has been shown that it is possible to discriminate different quantum computers by looking at the (unknown) noise fingerprints that characterize each device. Thus, what the ML techniques—presented here—could provide as an added value is to witness whether such noise fingerprints are time-correlated, and possibly how much the time-correlation is non-Markovian. For such experimental noise benchmarking, as in [61], it could be convenient to fix the connections among quantum gates (i.e., the underlying topology), and then consider more realizations of the implemented quantum dynamics affected by noise, so as to avoid the monitoring of the dynamics (cf. section 3). Notice that to make a successful classification among i.i.d. and (non-)Markovian noise samples, one should be able to previously label them as Markovian or not Markovian, and more in general to understand the main features of the noise acting on the machines. However, such task is usually very hard to be carried out. Thus, as a more plausible strategy, we propose to compare the experimental data to the theoretical prediction at the level of multi-time measurement statistics, according to the following two steps:

• (1)  
  Discriminate if and how much the measurement statistics (provided by the distributions ${{ \mathcal P }}_{{t}_{k}}$ with k=1, ... ,M, which have been measured on the real quantum devices on multiple runs) differ from the corresponding theoretical predictions. Such a difference, here on denoted as ${\mathfrak{D}}$, between theoretical and experimental data returns an effective prediction of the presence of noise on the machines.
• (2)  
  Conditionally to step 1), evaluate with ML-models the presence or the absence of functional relations ${ \mathcal F }$ that link the difference distributions ${{\mathfrak{D}}}_{{t}_{k}}$ in correspondence of the time instants t_k . If two consecutive instances of ${\mathfrak{D}}$ at times t_k−1 and t_k are no functionally related, then the noise is originated from a i.i.d. stochastic process. If, instead, there exist a functional ${{ \mathcal F }}_{{t}_{k-1}:{t}_{k}}=f({{\mathfrak{D}}}_{{t}_{k-1}},{{\mathfrak{D}}}_{{t}_{k}})$ that links together ${{\mathfrak{D}}}_{{t}_{k-1}}$ and ${{\mathfrak{D}}}_{{t}_{k}}$ in a non-trivial way, then the noise would come from a Markovian process. Finally, the noise process would be non-Markovian for functional relations ${{ \mathcal F }}_{{t}_{k-n}:{t}_{k}}$ defined over multi times with n > 1.The empirical characterization of ${ \mathcal F }$, as well as the witness of non-Markovianity in the noise samples, can be obtained through the use of generative models. We can train three different models: (i) one model that generates ${{\mathfrak{D}}}_{{t}_{k}}$ by processing ${{\mathfrak{D}}}_{{t}_{k-n}},\,\ldots \,,{{\mathfrak{D}}}_{{t}_{k-1}};$ (ii) another generative model that returns ${{\mathfrak{D}}}_{{t}_{k}}$ by taking ${{\mathfrak{D}}}_{{t}_{k-1}}$ as input; (iii) a model that directly generates ${{\mathfrak{D}}}_{{t}_{k}}$. If these three models have the same accuracy, then the noise process is likely i.i.d.. While, if the generative model (ii) has higher accuracy with respect to (iii) and the same accuracy than (i), then the noise process would be Markovian. Finally, if (i) has higher accuracy with respect to the models (ii) and (iii), then the noise process is non-Markovian.

To conclude, according to our proposal that will be tested in a forthcoming paper, time-correlations in the noisy samples of the distributions ${{ \mathcal P }}_{{t}_{k}}$, with k=1, ... ,M, can be determined by classifying functional relations ${ \mathcal F }$ linking the difference distributions ${\mathfrak{D}}$, obtained by comparing the theoretical and measured values of ${ \mathcal P }$ for a set of time instants. This is equivalent to discriminate coloured noise processes originated by different discrete Markov chain with non-zero transition matrices T. As a remark, it is still worth noting that also the experimental realization of the proposed procedure can be performed on multiple runs without the need to implement sequential measurements routines. Hence, each time a projective measurement is performed and the resulting measurement outcome recorded, the implemented (noisy) quantum circuit shall be executed from the beginning.

As outlook, we plan to test the ML-models employed in this paper on reconfigurable experimental platforms as the ones in [63, 64], even affected by multiple noise sources. Moreover, we also aim to adapt our ML methods (and especially ANNs) to reconstruct noise processes with time-correlation as key feature in the context of regression task instead of classification. Indeed, our proposal is to provide accurate estimates of both the probability distribution Prob(g) and the transition matrix T, and the analysis would be extended for the prediction of spatially-correlated noise sources. In this way, ML approaches would represent a very promising, and possibly more accurate, alternatives to other noise-sensing techniques, e.g., those recently discussed in [65, 66].

A well-known problem in ML is the generalization to data shift. A model that is trained on a data set sampled from a specific data distribution will work correctly only with data sampled from the same distribution. In this paper, we used only synthetic data to evaluate the correctness of the training process and the ML techniques. Thus, in order to validate this approach to real data, we should first collect them. This is out of the scope of the current work, but, as a remark, we can delineate three possible ways to build a real experimental data set. The first strategy is to acquire information a-priori on noise sources affecting the quantum system of interest in some experimental contexts by means of standard spectroscopy techniques, so that we can train the proposed ML-models to discriminate between unseen classes of noise. In this way, the initial effort in building a training data set that also contains experimental data is counterbalanced by the possibility to predict noise features by means of faster classification tasks. The second strategy, which has been employed in [61], is to collect experimental data that comes from different noisy measurement statistics whose noise processes are not necessarily known. The ML-models, then, are trained to classify (unknown) noise sources in distinct unseen sets of measurements. Finally, the third strategy, which is aimed to reduce the effort in building an informative experimental data set, is to train the ML-models first on synthetic data and then to fine tune the training on a smaller further data set with only experimental data. In such a case, it is beneficial to adopt a synthetic data set that closely adapts to the real experimental setup. For instance, a simulated extra error can be added to the measurement statistics (in our paper provided by the distributions ${{ \mathcal P }}_{{t}_{k}}$, k = 1, ... ,M) to take into account the finite number of measurement shots used for their computation.

Support Vector Machine (SVM) [33] is one of the first ML-model originally used to carry out classification tasks. The SVM training consists in finding the hyperplane that separates the elements x_q in two groups: one with the label y_q = 1 and the other with y_q = −1. The final hyperplane, solution of the classification, is the one having the maximum geometrical distance from the two parallel hyperplanes that are defined by the subsets of x_q called the support sets. When the data is not linearly separable, the kernel trick allows to increase the dimension of the features space in a way that the data becomes linearly separable in the new space.

Historically, SVM is a generalisation of the Support Vector Classifier (SVC) that, in turn, is an improved version of the Maximal Margin Classifier (MMC) [33]. MMCs aim at finding the hyperplane separating the two aforementioned classes of points, such that the distance between the hyperplane and the nearest points of the classes (commonly denoted as margin) is maximised. If the points of the data set are not linearly separable, then the value of the margin is negative. In such a case, the MMCs cannot be adopted. SVCs increase the performance of MMCs, by allowing some points of the data set, called slack variables, to be in the opposite part of the hyperplane with respect to the others of the belonging class. If the data set exhibits a non-linear bound between the two classes of points, SVCs are not able to correctly separate them, albeit the method returns a solution. Finally, SVMs extend the capabilities of SVCs by increasing the number of dimensions of the feature-space, such that in the new space the data set becomes linearly separable.

There are several classification problems (as for example the ImageNet Large Scale Visual Recognition Challenge [56] employing millions of images with hundreds of categories) that are solved through SVM but with a quite high residual classification error. For this reason, in order to improve the performance in solving classification problems, ANNs have been recently (re-)introduced as more-performing tools, and since 2012 have been extensively used [47, 56–58].

A MLP is composed of a variable number of fully connected layers, each of them with a variable number of artificial neurons. A single artificial neuron with I inputs (x) calculates the output as

$$\begin{eqnarray*}&&\hat{y}\equiv \sigma ({{\bf{w}}}^{T}\cdot {\bf{x}}+b)\end{eqnarray*}$$

that is the weighted sum of the inputs ${\bf{x}}\in {{\mathbb{R}}}^{I}$ with weights ${\bf{w}}\in {{\mathbb{R}}}^{I}$, plus a bias term $b\in {\mathbb{R}}$, followed by a nonlinear activation function $\sigma :{\mathbb{R}}\to {\mathbb{R}}$. The most common activation functions σ( · ) are: The sigmoid $\sigma (x)\equiv {\left(1+{e}^{-x}\right)}^{-1};$ the hyperbolic tangent $\sigma (x)\equiv \tanh (x);$ and the rectifier $\sigma (x)\equiv \max (0,x)$ [67, 68]. A single MLP-layer, composed of O neurons with I inputs, calculates

$$\begin{eqnarray*}&&\hat{{\bf{y}}}\equiv \sigma ({W}^{T}\cdot {\bf{x}}+{\bf{b}}),\end{eqnarray*}$$

where $\hat{{\bf{y}}}\in {{\mathbb{R}}}^{O}$ is the output vector, $W\in {{\mathbb{R}}}^{I\times O}$ is a matrix that collects all the weight vectors of the single neurons, and ${\bf{b}}\in {{\mathbb{R}}}^{O}$ is the vector of the biases. Finally, an MLP with L layers calculates

$$\begin{eqnarray}&&{\bf{h}}[l]\equiv \sigma \left(W{\left[l\right]}^{T}\cdot {\bf{h}}[l-1]+{\bf{b}}[l]\right)\end{eqnarray} \tag{ 6 }$$

with l = 1, ... ,L (index over the layers) and h[0] ≡ x. Thus, $\hat{{\bf{y}}}\equiv {\bf{h}}[L]$ is the output of the MLP, where W[l] and b[l] are, respectively, the weights and the biases of the l-th layer. Also the activation function may change depending on the specific layer. More concisely, the MLP can be denoted by the function

$$\begin{eqnarray}&&\hat{{\bf{y}}}=f({\bf{x}};\theta ,\xi )\end{eqnarray} \tag{ 7 }$$

of the inputs x. The function f is parametrised by the set θ ≡ {W[1], b[1], ... , W[L], b[L]} and by the fixed hyperparameters ξ defining the number, the dimension, and the activation functions of the MLP layers.

Let us now introduce the supervised learning process. For the sake of clarity, we just refer to the training of the MLP; however, the same notions can be applied in general to the supervised learning of vast majority of ANNs.

Equation (7) behaves like a generic function approximator [69]. Ideally, in the training process we would like to find the parameters

$$\begin{eqnarray}&&{\theta }^{* }=\mathop{\arg \,\min }\limits_{\theta }{L}_{{ \mathcal D }}(\theta ,\xi )\end{eqnarray} \tag{ 8 }$$

that minimise the theoretical risk function

$$\begin{eqnarray}&&{L}_{{ \mathcal D }}(\theta ,\xi )\equiv {{\mathbb{E}}}_{({\bf{x}},{\bf{y}})\sim { \mathcal D }}\left[{\ell }\left(f({\bf{x}};\theta ,\xi ),{\bf{y}}\right)\right],\end{eqnarray} \tag{ 9 }$$

i.e., the expected value of ℓ for (x, y) sampled from the distribution ${ \mathcal D }$ that generates the data set [32]. In equation (9), ${\ell }:{{\mathbb{R}}}^{O\times O}\to {{\mathbb{R}}}^{+}$ denotes the loss function (usually taken as a differentiable function, apart removable discontinuities) that measures the distance between the prediction $\hat{{\bf{y}}}$ and the desired output y. In general, the distribution ${ \mathcal D }$ is unknown; thus, the minimisation problem in equation (8) cannot be neither calculated nor solved. Indeed, one can dispose of a finite set $S=\{{\left({\bf{x}},{\bf{y}}\right)}_{1},\,\ldots \,,{\left({\bf{x}},{\bf{y}}\right)}_{n}\}$ of samples, to train, validate and test the ML-model. By considering the partition $\{{S}_{{tr}},{S}_{{va}},{S}_{{te}}\}$ of S, the theoretical risk function is approximated by the empirical risk function

$$\begin{eqnarray}&&{L}_{{S}_{{tr}}}(\theta ,\xi )=\displaystyle \frac{1}{| {S}_{{tr}}| }\displaystyle \sum _{({\bf{x}},{\bf{y}})\in {S}_{{tr}}}{\ell }\left(f({\bf{x}};\theta ,\xi ),{\bf{y}}\right)\end{eqnarray} \tag{ 10 }$$

that is the arithmetic mean of the loss function ℓevaluated on all the samples of the training set S_tr [32]. By minimising the empirical risk function ${L}_{{S}_{{tr}}}(\theta ,\xi )$ with respect to θ, the MLP is trained and θ^* obtained. Then, the validation set S_va is used to compute the empirical risk ${L}_{{S}_{{va}}}({\theta }^{* },\xi )$ that takes as input the optimal parameters attained by the minimisation of ${L}_{{S}_{{tr}}}$ (training stage). This procedure allows to check if the ML-model works also for unseen data. Notice that the minimisation of the training risk function ${L}_{{S}_{{tr}}}(\theta ,\xi )$ with respect to θ is performed step-by-step over time. After each step (also called epoch), the validation risk ${L}_{{S}_{{va}}}({\theta }^{* },\xi )$ is evaluated, and the minimisation procedure is stopped when the time-derivative of ${L}_{{S}_{{va}}}({\theta }^{* },\xi )$ becomes positive for several epochs, thus showing overfitting [70]. In case such time-derivative remains negative or constant over time, the procedure is ended after a predefined number of epochs. The validation set S_va can be also used to explore other configurations ξ of the ML-model: this process is called hyperparameters optimization. In particular, after completing the training procedure using two different set of hyperparameters ξ and $\xi ^{\prime}$, we obtain two minima θ^* and $\theta {{\prime} }^{* }$, and then compare ${L}_{{S}_{{va}}}({\theta }^{* },\xi )$ with ${L}_{{S}_{{va}}}(\theta {{\prime} }^{* },\xi ^{\prime} )$ to also choose the best hyperparameter. Finally, we use the test set S_te to calculate a significant metric (in our case, the classification accuracy) and report the results.

Regarding the hyperparameters optimization, it can be performed in different ways. The most basic technique is called grid search whereby the training and validation are carried out on a specific set of hyperparameters configurations. The random grid search considers configurations where each hyperparameter is randomly chosen within an a-priori fixed range of values. It has been proved to be more efficient than standard grid search [71]. A more sophisticated class of optimization methods is the Bayesian optimization [72] that updates, after the training of each hyperparameters configuration, a Bayesian model of the validation error. The best hyperparameters configuration is thus chosen as the one allowing for the lower guess validation error.

The most used optimisation algorithm to minimise equation (10) is the Stochastic Gradient Descent (SGD) [73–75] and its adaptive variants, such as Adaptive Moment Estimation (ADAM) [76], that changes the value of the learning rate η (i.e., the descent step) at each iteration. After having calculated the predictions $\hat{{\bf{y}}}$, the loss function ${\ell }(\hat{{\bf{y}}},{\bf{y}})$ is propagated backwards (backpropagation) in the ANN and its gradient in the weight space is calculated. Overall, the optimisation process consists in iteratively updating the value of the weights θ according to the relation

$$\begin{eqnarray*}&&{\theta }_{i}={\theta }_{i-1}-\eta {{\rm{\nabla }}}_{\theta }{L}_{{S}_{b}}({\theta }_{i-1},\xi ),\end{eqnarray*}$$

where i is the index for the descent step and ${S}_{b}\subset {S}_{{tr}}$ denotes the b-th set of samples, taken from the training set and used for the computation of the gradient. If ${S}_{b}\equiv {S}_{{tr}}$, the algorithm is called batch SGD; if S_b contains only one element is called on-line SGD; finally, the most common approach (we use it here) is mini-batch SGD that consider ∣S_b ∣ = B with B a fixed dimension [75]. Hence, the update of θ follows the descent direction of the gradient, with a magnitude determined by the learning rate η.

Now, let us introduce the specific loss function ℓconsidered in this paper. For classification problems with two or more classes, a common choice for ℓis the categorical cross entropy, which is defined as

$$\begin{eqnarray}&&{\ell }(\hat{{\bf{y}}},{\bf{y}})=-\displaystyle \sum _{j=1}^{O}{y}^{(j)}\mathrm{log}{\hat{y}}^{(j)}.\end{eqnarray} \tag{ 11 }$$

This function measures the dissimilarity between two or more probability distributions. Thus, to properly use the categorical cross entropy, it is convenient to choose the desired outputs y as Kronecker delta functions centered around the indices associated with each class to be classified. The model output $\hat{{\bf{y}}}$, instead, is normalised so that it represents a discrete probability distribution, i.e., a vector of positive elements summing to 1. This operation is obtained by using softmax [46] as the activation function of the last layer:

$$\begin{eqnarray}&&{\sigma }^{(i)}({\bf{z}})\equiv \displaystyle \frac{{e}^{{z}^{(i)}}}{{\sum }_{j=1}^{O}{e}^{{z}^{(j)}}}\end{eqnarray} \tag{ 12 }$$

where σ(z) is the vector having as elements σ⁽ⁱ⁾(z), with i = 1...,O, and z denotes the output of the last layer before the activation function.

In the experiments, the activation functions for the hidden layers of the MLP have been chosen among the sigmoid, hyperbolic tangent and rectifier functions accordingly to the hyperparameters optimization.

A Recurrent Neural Network (RNN) is an ANN specialised for sequence processing when the data set is expressed as

$$\begin{eqnarray}&&S=\left\{{\left(\left\{{{\bf{x}}}_{1},\,\ldots \,,{{\bf{x}}}_{{\tau }_{1}}\right\},{\bf{y}}\right)}_{1},\,\ldots \,,{\left(\left\{{{\bf{x}}}_{1},\,\ldots \,,{{\bf{x}}}_{{\tau }_{n}}\right\},{\bf{y}}\right)}_{n}\right\},\end{eqnarray} \tag{ 13 }$$

where τ_r defines the number of elements of the r-th sequence. RNNs can be used in tasks regarding Natural Language Processing (NLP) [48, 77–79], time series analysis [80] and, in general, all the tasks involving ordered set of data [81]. Note that, in general, for sequence-to-sequence problems also y can be a sequence of elements, as for example in machine translation where the inputs and outputs of the RNN are sentences in different languages [82]. In this paper, sequence-to-sequence problems are not considered, and we thus consider τ₁ ≡ ... ≡ τ_n ≡ τ.

A RNN is defined by the recurrent relation

$$\begin{eqnarray}&&{{\bf{h}}}_{t}=r\left({{\bf{x}}}_{t},{{\bf{h}}}_{t-1};{\theta }_{r},\xi \right)\end{eqnarray} \tag{ 14 }$$

where t ∈ {1, ... , τ}, ${{\bf{h}}}_{t}\in {{\mathbb{R}}}^{d}$ is a d-dimensional vector with d being an hyperparameter belonging to ξ and h₀ = 0 (vector of zeros). The recurrent relation (14) defines τ hidden representations h_t (to be seen as a memory) of the input sequence ${\{{{\bf{x}}}_{1},\,\ldots \,,{{\bf{x}}}_{\tau }\}}_{q}$ with q = 1, ... ,n. If the function r is implemented as an MLP (7) that takes as input the concatenation x_t ⊕ h_t−1 (usually called 'vanilla RNN'), the model suffers the so-called vanishing gradient problem [83, 84] such that the weights of the last layers of the RNN are updated only with respect to the more recent input data. The vanishing gradient problem occurs when the backpropagation is performed on a high number of layers, as it could happen in our case with a large value of τ (thus meaning long input sequences). In this regard, to mitigate the vanishing gradient problem, LSTM [85] and GRU [86] have been introduced. These methods use learned gated mechanisms, based on current input data and previous hidden representations, to control how to update the current hidden representation h_t . Specifically, if LSTM is used, equation (14) needs to be slightly modified as

$$\begin{eqnarray}&&{{\bf{s}}}_{t}=v\left({{\bf{x}}}_{t},{{\bf{s}}}_{t-1};{\theta }_{v},\xi \right)\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{\bf{h}}}_{t}=r\left({{\bf{x}}}_{t},{{\bf{s}}}_{t},{{\bf{h}}}_{t-1};{\theta }_{r},\xi \right)\end{eqnarray} \tag{ 16 }$$

where s₀ = 0 and v, r are, as usual, nonlinear functions. Both for GRU and LSTM, the nonlinearity of the recurrent relations is due to the adoption of the hyperbolic tangent and sigmoid functions, where the latter are employed only for the gating mechanism. It is worth noting that in equation (15) s_t is a state vector that allows to differently propagate over time specific elements of the hidden representations h_t depending on the input data. This means that, at any time t, the hidden representation h_t depends not only on the input x_t and the previous hidden representation h_t−1 but also on the state vector s_t . For further details, refer to [48, 85, 86].

RNNs are usually considered deep-learning models, due to the high number of layers, when they are unfolded on the sequence dimension for t = 1, ... , τ. The key aspect of deep learning is the automatic extraction of features by means of the composition of a large number of layers; an increasing (deep) number of layers is typically used to extract features with increasing complexity [87].

Moreover, RNNs can be extended considering more layers [88] and processing the data bidirectionally [89]. Regarding the latter, one can define two sets of hidden representations: One for the forward and the other for the backward direction, where the t-th hidden representation depends respectively on the (t − 1)-th or (t + 1)-th one. More formally,

$$\begin{eqnarray}&&{{\bf{h}}}_{t}[l]=r\left({{\bf{h}}}_{t}[l-1],{{\bf{h}}}_{t-1}[l];{\theta }_{r}[l],\xi \right)\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\widetilde{{\bf{h}}}}_{t}[l]=r\left({\widetilde{{\bf{h}}}}_{t}[l-1],{\widetilde{{\bf{h}}}}_{t+1}[l];{\theta }_{\tilde{r}}[l],\xi \right)\end{eqnarray} \tag{ 18 }$$

with l=1, ... ,L and ${{\bf{h}}}_{t}[0]\equiv {\widetilde{{\bf{h}}}}_{t}[0]\equiv {{\bf{x}}}_{t}$.

Now, let us explain how to use the hidden representations to calculate the prediction $\hat{{\bf{y}}}$ in output from the ML-model. The common approach to calculate the prediction in classification problems is to use the RNN as an encoder of the sequence and to scale the dimension of the last hidden representation ${{\bf{h}}}_{\tau }[L]\oplus {\widetilde{{\bf{h}}}}_{1}[L]$ (in the more general case of bidirectional models) to the one of the output vector. This scaling can be done through a fully connected layer, or, more in general, by means of an MLP, i.e.,

$$\begin{eqnarray}&&{\bf{a}}\equiv {{\bf{h}}}_{\tau }[L]\oplus {\widetilde{{\bf{h}}}}_{1}[L]\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\hat{{\bf{y}}}=f\left({\bf{a}};{\theta }_{f},\xi \right).\end{eqnarray} \tag{ 20 }$$

Then, we can use SGD to minimise an empirical risk function similar to equation (10) of MLPs.

It is possible to consider different forms of aggregation a for the hidden representations h_t [L], with t = 1, ... , τ, instead of using only the last hidden representation as in equation (19). In this regard, attention mechanisms [90–92], also in hierarchical forms [93], perform a weighted average of the h_t [L] where the weights are learned together with the ML-model. In detail, equation (19) becomes:

$$\begin{eqnarray}&&{{\bf{u}}}_{t}={{\bf{h}}}_{t}[L]\oplus {\widetilde{{\bf{h}}}}_{t}[L]\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{\bf{v}}}_{t}=\tanh \left({{\bf{W}}}^{T}\cdot {{\bf{u}}}_{t}+{\bf{b}}\right)\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\alpha }_{t}\equiv \displaystyle \frac{{e}^{\langle {{\bf{v}}}_{t},{\bf{c}}\rangle }}{{\sum }_{j=1}^{\tau }{e}^{\langle {{\bf{v}}}_{j},{\bf{c}}\rangle }}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{\bf{a}}\equiv \displaystyle \sum _{t=1}^{\tau }{\alpha }_{t}{{\bf{u}}}_{t},\end{eqnarray} \tag{ 24 }$$

where 〈·, · 〉 denotes the dot product and c is a learned vector that is randomly initialised and jointly learned during the training process as in [90–93]. Another form of aggregation a is the max pooling aggregation, whereby each element a^(j) of a just refers to a single value of t. In this case, equation (19) equals to

$$\begin{eqnarray}&&{{\bf{a}}}^{(j)}=\mathop{\max }\limits_{t}\,{{\bf{u}}}_{t}^{(j)}.\end{eqnarray} \tag{ 25 }$$

where the expression of u_t is provided by equation (21). In this way, each element ${{\bf{u}}}_{t}^{(j)}$ of the hidden representations (for t = 1, ... , τ) learns to detect specific features of the input data within all the interval [1, τ].

Finally, another approach, which we do not use here, is to consider the RNN as a transducer that produces an output sequence ${\hat{{\bf{y}}}}_{t}$ for $t=1,\,\ldots \,,\widetilde{\tau }$ (generally $\widetilde{\tau }\ne \tau$) in correspondence of the input sequence x_t with t = 1, ... , τ [48, 82, 94].

All the ML-models are realized in PyTorch and have been trained on the six different data sets using a DELL® Precision Tower workstation with one NVIDIA® TITAN RTX® GPU with 10 Gb of memory, 88 cores Intel® Xeon® CPU E5-2699 v4 at 2.20 GHz and 94 Gb of RAM.

We train the ANN models in mini-batches of dimension 16 by means of the SGD using ADAM [76] and learning rate η = 10⁻³. We optimize the hyperparameters ξ with ASHA [95] as scheduler and Hyperopt [96, 97] (Hyperopt belongs to the family of Bayesian optimization algorithms) as search algorithm in the framework Ray Tune [98]. For the MLP models, the hyperparameters optimization defines: (i) the activation functions to be used, (ii) the number of layers, and (iii) their dimension, within the following search space: $\sigma \in \{\mathrm{relu},\mathrm{sigmoid},\tanh \}$, L ∈ {2, 3, 4, 5, 6} and $\dim ({\bf{h}}[1])\equiv ...\equiv \dim ({\bf{h}}[L])\in \{d\in {\mathbb{N}}\,| \,1\leqslant d\leqslant 512\}$. Instead, for the RNN models the search space is L ∈ {1, 2, 3, 4} for the number of recurrent layers (L ∈ {1, 2, 3, 4, 5, 6} for the NM task with t₁₅ = 0.1), and $\dim ({{\bf{h}}}_{t}[1])\equiv \dim ({{\bf{h}}}_{t}[1])\equiv ...\equiv \dim ({\widetilde{{\bf{h}}}}_{t}[L])\equiv \dim ({\widetilde{{\bf{h}}}}_{t}[L])\,\in \{d\in {\mathbb{N}}\,| \,1\leqslant d\leqslant 512\}$ for the layers dimension. Regarding the ML-models m-biGRU-att andm-biLSTM-att, the search space includes also the dimension of the attention layer as in equations (22) and (23), i.e., $\dim ({\bf{c}})\equiv \dim ({{\bf{v}}}_{1})\equiv ...\equiv \dim ({{\bf{v}}}_{\tau })\in \{d\in {\mathbb{N}}\,| \,1\leqslant d\leqslant 512\}$. In the hyperparameters optimization of all the MLP and the RNN models, we have also used regularization methods as weight decay [99] and dropout [100]. They are able to mitigate overfitting; in particular, the former adds a penalty (chosen among {0, 10⁻⁴, 10⁻³}) to the risk function L_S (θ, ξ) with the aim to discourage large weights. Instead, using dropout, the outputs of the artificial neurons during the training are forced to zero with a probability among {0, 0.2, 0.5}.