Optimization of the Memory Reset Rate of a Quantum Echo-State Network for Time Sequential Tasks

Quantum reservoir computing is a class of quantum machine learning algorithms involving a reservoir of an echo state network based on a register of qubits, but the dependence of its memory capacity on the hyperparameters is still rather unclear. In order to maximize its accuracy in time--series predictive tasks, we investigate the relation between the memory of the network and the reset rate of the evolution of the quantum reservoir. We benchmark the network performance by three non--linear maps with fading memory on IBM quantum hardware. The memory capacity of the quantum reservoir is maximized for central values of the memory reset rate in the interval [0,1]. As expected, the memory capacity increases approximately linearly with the number of qubits. After optimization of the memory reset rate, the mean squared errors of the predicted outputs in the tasks may decrease by a factor ~1/5 with respect to previous implementations.


Introduction
The continuous development of quantum computers [1][2] [22] [24][23] [25] has found in machine learning one of its natural applications [3] [4]. In particular, the field of quantum machine learning (QML) has been widely explored mainly through hybrid quantum-classical algorithms which take advantage of both quantum and classical parts during the training. Many of such applications involve variational algorithms to solve classification problems [12] [11], data reconstruction [3] and generative adversarial networks [5]. Differently, quantum reservoir computing (QRC) deals primarily with sequential temporal tasks where the fading memory of the system plays a crucial role. Based on the concept of reservoir computing [7] originally embodied by echo state networks [21] and liquid state machines [9] respectively, quantum reservoir computing was first proposed in 2017 by Fujii and Nakajima [6]. There, an encoding of a quantum reservoir on a register of qubits is proposed. The implementation consisted of the simulation of a NMR machine [14,10,18]. Other proposals based on gate-model hardware have followed [15,13]. In QRC, instead of usual rate-neurons, a register of qubits provides the reservoir of an echo-state network [19,16,21,20], thus realizing a quantum ESN (qESN). By a suitable quantum evolution of the system, the reservoir expands non-linearly the input values in a high dimensional system so that it is sufficient to add a final linear readout to make the network capable of emulating any nonlinear map. Furthermore, the reservoir retains information about past inputs in its quantum state enabling the network to tackle temporal tasks for which memory is essential.
Quantum reservoir computing, and in particular quantum echo-state networks, appear as promising implementations of quantum machine learning for temporal tasks, but a detailed analysis of the dependence of the memory properties of a qESN on its hyperparameters is still largely unadressed. With the spirit of filling such gap, we study the memory dependence on the reset rate which regulates the probability with which the quantum reservoir resets its internal state during the evolution.
Previously, some of the authors explored QML through quantum annealers for unsupervised learning [3], quantum generative adversarial networks [5,8], quantum feed-forward neural networks [11] and quantum variational tensor networks for multiclass supervised learning [12]. Here we explore supervised learning of time-series by a qESN through an embodiment on gate model quantum computers. The implementation of the qESN used in our work is particularly suited for noisy intermediate-scale quantum (NISQ) computers. In particular we take in consideration the subclass of universal QRs which evolves with a resettable convex linear combination of complete trace preserving (CPTP) maps, introduced in Ref. [15].
We investigate the memory of the reservoir by employing the memory capacity (MC) metric [6]. The memory capacity quantifies the short-term memory of the reservoir by testing its ability to reproduce previously received input values by using its current internal state. In particular, in Ref. [6] an empirical evaluation suggested that a qESN of 5 qubits could exhibit a higher memory capacity than a classical ESN of 500 nodes.
In our experiments we measured the memory capacity for different values of the reset rate , a fundamental hyperparameter of the QR subclass studied in this work, which regulates the rate at which the quantum state of the reservoir is reset during the time evolution, thus driving its fading memory behaviour. We compare the memory capacity of different reservoirs based on 3, 5 and 7 qubits, respectively.
Furthermore, the value of influences also the overall performance of the network. We tested our qESN over three different tasks, all of them requiring the emulation of a nonlinear map with fading memory, that is with outputs which depend increasingly less on previous inputs as time flows. We employed two, relatively simple, second-order maps and a more involved realistic map consisting in a pair of differential equations related to the simplified dynamics of an aircraft. For all of the tasks we injected into the system random input variables and trained the network to predict the exact output values for input data not used during the training, according to the specific map. For each task we tested the implementation on a simulator with different values of . Next we selected the optimal one to run the experiments on quantum hardware. The qESN is realized via quantum circuits on both the Qiskit simulator and the IBM ibmq_casablanca quantum hardware.
The results of the memory investigation show that the MC is maximized for the central values of ∈ [0, 1], while it sensibly decreases for > 0.5. Furthermore for the optimal = 0.5 the MC increases approximately linearly with the number of qubits of the reservoir, confirming the result obtained in Ref. [6] with a different evolution for the reservoir by using a simulated NMR machine.
The qESN is able to predict outputs which match those expected with low errors, measured as normalized mean squared errors (NMSE), for all the three maps on which it was tested. In particular, by choosing the optimal of the network for each task, lower errors are obtained in comparison with previous implementations in Ref. [15].
The lowest error in the tasks is obtained in correspondence of the values for which the MC is maximum.
The paper is organized as follows. In Section 2 we summarize the theory behind QRC and in particular quantum echo-state network while in Section 3 we describe the actual implementation as a circuit on a gate model quantum computer. The results are shown in Section 4 along with the detailed description of the tasks and the comparison between different settings for the reservoir with various values of the reset rate and numbers of qubits. The conclusions are drawn in the last Section 5.

Method
The key idea of QRC resides in using a N − qubits quantum register as reservoir for an ESN. In such a way it is possible to exploit the high dimensionality of the Hilbert space associated to N qubits, which effectively acts as a reservoir with a large number of nodes. Specifically, the nodes of the reservoir can be associated to the basis elements of the operator space, which for a N qubits system amount to 4 N elements. We consider as basis elements the N -qubit Pauli operators defined as: where {σ 01 , σ 10 , σ 11 } = {X, Y, Z} are the three Pauli matrices while σ 00 = I, the 2 × 2 identity matrix, and i is a natural number described by a binary string ranging from 000...0 to 111...1, build from the indexes of the Pauli matrices in the tensor product. For a N qubit reservoir, a single operator P i is the tensor product between N matrices σ ∈ {σ 01 , σ 10 , σ 11 , σ 00 }. There are 4 N possible total combinations of such matrices, so that i = 1, 2, ...4 N . Thanks to the orthogonality property a generic density matrix ρ describing the quantum state of the system can be written as where the real numbers can be regarded as classical coordinates, or nodes, of our qESN. Assuming P 4 N = N k=1 I k , one observes that the last element of Eq. (4) is constrained by normalization, i.e. x 4 N = 1. Therefore, for a system of N qubits there are 4 N − 1 internal nodes of the reservoir.
We are interested in the learning of a temporal tasks where the input signal is an ordered array {u(k)} L k=1 of L real values with u(k) ∈ [0, 1] for every k. The input can be injected into the quantum reservoir with a pure states encoding, for example by initializing at each time step the first qubit of the reservoir in the state |ψ 1 = u(k)|0 + 1 − u(k)|1 [6]. Instead, we follow the approach of Ref. [15] so to encode the input as a classical mixture. Such an approach corresponds to initializing an ancilla qubit in the diagonally mixed state: The ancilla qubit is then coupled with the N −qubits register of the reservoir in order to evolve the system with an input dependent map T acting on the density matrix ρ R of the reservoir, initialized at t = 0 in the pure state ρ R (0) = |0 0| ⊗N . In particular, at each time step t = k the density matrix ρ R evolves as: where T (u(k)) is a map which acts within the space of density matrices (i.e. T : C N × C N → C N × C N ) whose structure is described in Section 3. The time evolution of the reservoir serves the purpose of nonlinearly expanding the input signal, which is a sequence of scalar values, in a higher dimensional space, i.e. the reservoir degree. Indeed, by choosing a suitable map T (u(k)), the output signals extracted from the reservoir are nonlinear functions of the input values [6], [14]. In this way the network is able to emulate a nonlinear map by just adding a linear readout at its output. After each injection of an input value and one-step evolution of the reservoir state, N signals are extracted from the system and passed to the linear readout layer. The expectation values of the Z operators i.e.
provide the output signals. Such averages are computed as Z i = Tr[Z i ρ] and correspond to N so-called true nodes of the reservoir. The remaining 4 N − 1 − N so-called hidden nodes are never measured. In Figure 1 the structure of the qESN employed in the experiments is shown. In the next subsection the process of training the qESN and how the weights are updated is explained.

Training
In the following, the goal of the temporal task is to learn a non linear Since in general the output at a given time depends also on the history of the input signal at previous times, that , some kind of memory mechanism is required on the network. Indeed the quantum evolution in Eq. (6) allows the system to retain information about past inputs, so it provides the ability to correctly emulate the given map f . In the following, three different maps are considered, all of them with fading memory, i.e. each output depends increasingly less on previous inputs as time proceeds. Further details on them are given in Section 3.1.
The basic idea of reservoir computing is to make a clear-cut distinction between the reservoir and the linear readout and train them differently. Specifically, once all hyperparameters are properly set, the only weights which are adjusted during the training in an ESN are those of the linear readout, through simple linear regression.
Here, for every k th time step, N signals are extracted from the reservoir: let's denote by x ki the value on the i−th true node at the time t = k, with i = 1, 2, ..., N and k = 1, ..., L, where L is the total temporal length of the input signal. The x ki form the L × N design matrix X. As target values for the training and testing we use the {ȳ(k)} L k=1 signal obtained from the maps for the corresponding input {u(k)} L k=1 . These target values can be collected in a L × 1 array Y. So, at the end, linear regression is performed over the matrix equation Y = XW calculating the N trained weights as: where X † = (X T X) −1 X T is the Moore-Penrose pseudo-inverse of X. Once the weights are trained, the prediction values for k > L are computed by: whereX collects the output values when k = L + 1, L + 2 . . .. We recall that these output signals are read from a reservoir whose internal dynamics is not subject to any training procedure, as appropriate to an Echo State Network.

Experiments
In this Section we discuss the details of the actual implementation of the qESN. For the quantum part of the network, i.e. the reservoir, we employed both a 5 and a 7 qubit system, the former in order to compare the results with those in Ref. [15]. We implemented it as a circuit on both the IBM Qiskit simulator and quantum hardware ibmq_casablanca. For the output layer, instead, a classical linear readout with N input and one output nodes was used, which we implemented by linear regression using the Python library sklearn.
As discussed in Section 2 at every time step the system evolves under a map T (u(k)) which depends on the input value u(k) at that time step t = k. For the sake of consistency with previous literature, the stochastic evolution map T proposed in Ref. [15], acting on the density matrix of the reservoir, is considered, by: where "w.p." is a short hand for "with probability". Here 0 ≤ u(k) ≤ 1 is the input at the time t = k while T 0 and T 1 are two arbitrary completely positive trace-preserving (CPTP) maps. Hence they can be implemented with two unitary operators U 0 and U 1 according to the usual scheme for time evolution of density matrices in Quantum Mechanics: In the third alternative of Eq. (9), the quantum state is reset to some fixed   density operator σ with probability where 0 ≤ ≤ 1. The hyperparameter , which can be called reset rate, allows the qESN to have a fading memory. If = 1 there is no evolution and no dependence whatsoever on the inputs {u(k)} L k=1 . At the opposite end, when = 0, the network undergoes a mixture of unitary evolution steps. The net evolution is not unitary, but it might retain too much information on the initial state, hindering the driving capability of the inputs.
A particular history of the stochastic process in Eq. (9) corresponds to a quantum circuit where we encode the input on an ancilla qubit in a mixed state, as in Eq. (5). Then we couple the ancilla with the reservoir through controlled unitary operators U j as shown in Figure 1(b). To implement the evolution map (9) in Qiskit, which does not allow for mixed-state qubit preparation, we employ an ensemble description [15]. Specifically we reproduce the mixed state encoding of the input by preparing an ensemble of quantum circuits all in the same initial pure state |ψ(0) = |0 ⊗N (with N = 5 or N = 7) and evolving them with operators U 0 or U 1 with the corresponding probability law. In particular we prepare an ensemble of N c circuits in the all spin up state and at every time step t = k we evolve each of them by applying either one of the two unitary operators U 0 or U 1 , or we reset the state |ψ(k) to |ψ(0) : The two unitary operators U 0 and U 1 can be chosen arbitrarily, provided they do not produce trivial dynamics, as would be for example if U 0 ∝ U 1 and do generate a nonlinear dependence on {u(k)} L k=1 . To ensure that the next results can be compared with prior art, the implementation follows the same choice already discussed in Ref. [15]. The circuits for the two unitary operators (shown in Figure 2) consist of U 0 as a highly non trivial operator composed only by CNOT gates, while U 1 takes into account 3D spin rotations on each qubit, respectively.
Summarizing, N c circuits of N −qubits (with N = 5 or N = 7) each are created and initialized, by default, at the state |0 ⊗N . At each time step t = k, for each one of such circuits, the sequence of gate U 0 or U 1 is applied accordingly to the input dependent probabilities as from Eq. (11). With a probability , instead, each circuit is reset in the initial state |0 ⊗N . For each one of the N c circuits, after every step evolution, the measurement of the Z j operators is done. Such process is repeated for all the L time steps and the results are averaged over all the N c circuits. We finally obtain the design L × N matrix X which collects all output histories from the reservoir. Such matrix is then used as argument of the linear readout to finally obtain the output signal of the ESN and compute the trained weights W trained with linear regression as in Eq. (7). In the next section, the details of the tasks used to test our qESN are described.

Description of the tasks
In this Section the three tasks over which the qESN has been tested are explained. The three tasks are taken from the ones employed in Ref. [15], in such way it was possible to compare the results. All of them consist in learning a nonlinear map, that means to induce the system to emulate a selected map predicting the correct outputs for input values not seen before.
For all the three tasks the input is a length L array of random values extracted in [0, 1] where the k −th entry is u(k) of from the previous Sections. In tasks I and II we employ two maps that process the input by updating linearly an internal highly dimensional state z at each time step k; then the map's output value is given by a polynomial readout h: Here z, c ∈ R n with n = 2000 and c an arbitrary fixed array, while h is a polynomial function, in our case taken to be of second degree, which returns a single scalar value. Specifically, c and all the parameters of the h function were selected randomly in [−1, 1], while the initial condition z(0) was chosen to be identically null. The two maps differ for the structure of the matrix A ∈ R n×n : while in the first case we used a dense matrix with a maximum singular value of σ(A) max = 0.5, in the second case a 95% sparse matrix with σ(A) max = 0.99 was employed. The third task instead relates to the simplified dynamics of an aircraft moving in the horizontal plane. In particular, if the velocity is constant then it can be shown [17] that such dynamics is defined by the couple of differential equation:ẋ Whereũ represents the rudder deflection at a given time which, as pointed out in Ref. [17], should be in the range [−0.5, 0.5]. So at each time t = k the relation betweenũ and the input values u(k) isũ = 0.5u(k) − 1. The two variables x 1 and x 2 , instead, concern the slideslip and yaw angle respectively. The output is then taken to beȳ(k) = x 1 (k) and the target values are obtained using the ode45 solver in MATLAB with a sampling time of t = 1. For all three tasks, the target signal is the output array with entries {ȳ(k)} L k=1 , given the input {u(k)} L k=1 . Specifically, we set L = 60 in the experiments on the Qiskit simulator while L = 30 on the quantum hardware for time taking reasons, as we will explain in the next Section.

Results
In this Section we show the results obtained for the three tasks both on the simulator and on the quantum hardware. Since quantum measurements perturb the state, after the extraction of the reservoir signal the information stored in the qubits is lost and the qubit is no longer usable. To circumvent this problem, after the measurement at t = k every circuit is reinitialized to |0 ⊗N and left to evolve unperturbed until t = k + 1, where a new measurement is performed. Therefore, to process an input signal of length L, one needs to reinitialize the circuit L − 1 times and for each time step to increase by 1 the number of applications of the evolution in Eq. (11), requiring a total of 1 + 2 + ... + L = L(L + 1)/2 applications of the gates in Eq. (11). Furthermore, every circuit at every time step is run for a total number S of shots. Since, as we said, for the Monte Carlo sampling we employ N c circuit copies, the total number of circuits runs and application of Eq. (11) are N c SL and N c S(L + 1)L/2 respectively. In order to reduce the computational time, we decided to take in consideration shorter input signal of L = 30 and fewer circuit copies N c = 300 for the experiments on the quantum hardware compared to the L = 60 and N c = 1024 employed on the simulator. All the parameters of the experiments are summarized in the Table 1.

Results on the simulator
For the experiments on the simulator we extracted L = 60 random values in [0, 1] which formed the entries of the input array. In particular, after discarding the first 10 values, the training was performed over the next 35 inputs, i.e. until the time step L T = 45. Such values were used to adjust the readout weights by linear regression. The last 15 values, which correspond to the time steps k = L T + 1, ..., L, were employed for the testing where we compared the output from the qESN to the {ȳ(k)} 60 k=46 expected outputs  obtained from the map of the given task. To asses the goodness of the performance we measured the discrepancy between the expectedȳ(k) and obtained y(k) output by computing a normalized mean square error (NMSE): where ∆ 2 y = L l=L T +1 (y(l) − µ) 2 with µ = 1/(L − L T ) L l=L T +1 y l . We performed each task for different values of ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}, always keeping the same input signal, in order to investigate which one performs better and then use it on the quantum hardware. The experiments were repeat 15 times with different inputs and the final errors of the testing were averaged. In Figure 3(a) we show the results obtained in a single experiment for each of the three tasks while in Figure 3(b) the averaged errors evaluated over the 15 different inputs are plotted. The lowest errors are obtained with = 0.5 for Task I and II while = 0.4 for Task III. Figure 3(b) shows the average errors for a reservoir made of 7 qubits. Its evolution operators U 0 and U 1 are based on the same structure of those used with the 5 qubits case and adapted them to a 7 qubits register, as shown in Figure 2. Except for Task III, the 7 qubits reservoir returns a lower error for every values. For all the three Tasks the optimal of the 7 qubits reservoir slightly differs from the 5 qubits case ( = 0.6 vs = 0.5 for Tasks I and II, = 0.5 vs = 0.4 for Task III). The following Section deals with quantum hardware with a reservoir of 5 qubits and therefore we selected the optimal value accordingly.

Results on quantum hardware
For the experiments on quantum hardware we employed the IBM quantum computer ibmq_casablanca. As shown in Table 1 we used an input of L = 30 values, in this case the training was performed over the first L T = 23 entries of the input array (the effective training was performed over 19 values as the first 4 inputs were discarded), while the last 7 were used for the testing. We run our experiment using an value of = 0.5 which is the best for Task I and II and still gives low error for Task III, as seen from Figure 3(b). The results are shown in Figure 4, while the performances of the network in each task are summarized in the Table 2. With respect to prior art, which is represented by the implementation of the same operators on the 5 qubits Ourense quantum computer of IBM at = 0.1 [15], we improve the MSE from 0.24 to 0.09 for Task I, from 0.68 to 0.12 for Task II, respectively, while the results on the Task III are comparable (MSE is 2.56 instead of 2.3), but one should keep in mind that we used less than one third of circuit copies.

Evaluation of the memory of the qESN
Finally, we present the dependence of the memory of the qESN on the reset rate . The memory of the system is quantitatively described by the memory capacity (MC) metric [6]. Such quantity determines the ability of the qESN to reproduce at a time t inputs received at the earlier time t − τ .
As τ increases the reservoir should gradually loose information about its past, so the accuracy of the predicted output should decrease. We run our  . The correlation decreases, as expected, as τ increases. Furthermore, at higher the C(τ ) decreases more rapidly. A slower decrease of C(τ ) can be seen for the 7 qubits reservoir. The histogram in (b) represents the MC at different for the 5 qubits reservoir (in blue). It is almost constant from = 0.2 to = 0.5, then it starts to decrease. The comparison with the 3 and 7 qubits reservoirs shows an increase of the memory capacity with the number of qubits. experiment on the Qiskit simulator. In particular we prepared an input array u = {u(k)} 60 k=15 with random values u(k) ∈ [0, 1] and we injected into the qESN the first 45 values as training set, picking as target array the same sequence shifted by a delay τ , i.e. u τ = {u(k − τ )} 60 k=15 . After discarding the first 14 output values, we trained the network by comparing output and target for the times k = 15, ..., 45, while the last 15 values were used for the testing as we did in the three tasks in Section 4.1. We then calculated the correlation between the output signal of the network y = {y(k)} 60 k=45 and the target t τ for a given τ : where COV(x, y) and σ(x) represent the covariance between x and y and the standard deviation of x, respectively. The memory capacity M c is then computed by summing the values of C(τ ) for different τ : where in our case τ max = 8. The results for the memory capacity tests are discussed in the next Section. In Figure 5(b) the histograms representing the memory capacity at different 's were displayed for reservoirs of 3, 5 and 7 qubits. In Figure 5(a), instead, the behavior of C(τ ) as τ varies is shown for various 's. To test the memory capacity, we employed three different reservoirs of 3, 5 and 7 qubits, respectively. For each one, the performance for different values of the reset rate were evaluated.
In the case of the 3 qubits reservoir we adopted for the evolution operator U 0 a circuit with one CNOT between the first-second and second-third quibits while for U 1 we took the same of Figure 2 considering only the first three qubits. In such way we were able to compare the memory of the three different reservoirs. The results in Figure 5 suggest that for the optimal = 0.5 the increase in the value of the MC (2.43, 3.09, 3.76, for N = 3, 5, 7 respectively) with the number of qubits is approximately linear, confirming previous results [6] with a different implementation of a quantum reservoir. Such result is expected as increasing the number of qubits will increase linearly the number of observables we measure in the output extraction process, i.e. the Z expectation values measured on each qubit.

Conclusions
We implemented a quantum Echo-State Network on a gate model quantum computer and tested it over three different tasks consisting to predict outputs of nonlinear maps, to maximize its performances. We investigated the memory of the reservoir by evaluating its memory capacity for different values of the reset rate . The memory capacity of the reservoir remains stable at its highest values for the central values of ∈ [0, 1], while it starts to decrease for > 0.5, as expected since the reservoir resets its internal state more often. After choosing the best parameter for the reservoir evolution based on the simulator at 0.5, the results obtained on the quantum hardware show low NMSE. In particular, by comparing the MSE values with those reported in literature at = 0.1, we observe an improvement of 62% and 82% in Task I and II respectively. For the task involving the pair of differential equations, a comparable error is found by using less than one third of circuit copies.For completeness, it should be noted that the results in Ref. [15] are obtained on a different machine (Ourense) than the one employed in this work (Casablanca) as Ourense is no longer available on the IBM system. The results obtained might, of course, be affected by the different nature of noise in the two machines. However, the decrease of MSE for central values of ∈ [0, 1] is not influenced by the nature of the machines and, so, by choosing an optimal better results than Ref. [15] are in any case expected. We confirm that also for the optimal = 0.5 the memory capacity grows approximately linearly with the number of qubits. Such result is consistent with that obtained in Ref. [6] where a different reservoir evolution was considered and the implementation was done by a simulated NMR machine. There is no direct proportionality between the memory capacity of the network and its errors in the tasks. However, the best performances are found for the values for which the memory capacity is highest while the errors in the tasks start to increase as the memory capacity decreases at the largest values of . Our investigation shows that, in order to maximize the performances of the qESN, it is convenient to tune the hyperparameters of the quantum reservoir such that the network develops a high memory capacity.
We also notice that the optimal found in our work should not in principle depend on the task in consideration. In fact the dynamics of the reservoir can be regarded as that of many chaotic oscillators with frequencies which spread over several orders of magnitude. In that sense, the optimal setting of the reservoir should not be overly influenced by the time dependence that the function to be emulated sets over input data. On the other hand, we do expect the weights of the readout layer do strongly depend on such time scale.