Density Matrix Emulation of Quantum Recurrent Neural Networks for
  Multivariate Time Series Prediction
Quantum Recurrent Neural Networks (QRNNs) are robust candidates to model and
predict future values in multivariate time series. However, the effective
implementation of some QRNN models is limited by the need of mid-circuit
measurements. Those increase the requirements for quantum hardware, which in
the current NISQ era does not allow reliable computations. Emulation arises as
the main near-term alternative to explore the potential of QRNNs, but existing
quantum emulators are not dedicated to circuits with multiple intermediate
measurements. In this context, we design a specific emulation method that
relies on density matrix formalism. The mathematical development is explicitly
provided as a compact formulation by using tensor notation. It allows us to
show how the present and past information from a time series is transmitted
through the circuit, and how to reduce the computational cost in every time
step of the emulated network. In addition, we derive the analytical gradient
and the Hessian of the network outputs with respect to its trainable
parameters, with an eye on gradient-based training and noisy outputs that would
appear when using real quantum processors. We finally test the presented
methods using a novel hardware-efficient ansatz and three diverse datasets that
include univariate and multivariate time series. Our results show how QRNNs can
make accurate predictions of future values by capturing non-trivial patterns of
input series with different complexities.
José Daniel Viqueira, Daniel Faílde, Mariamo M. Juane, Andrés Gómez, David Mera
Density Matrix Emulation of Quantum Recurrent Neural Networks for Multivariate Time Series Prediction

###### Abstract

Quantum Recurrent Neural Networks (QRNNs) are robust candidates to model and predict future values in multivariate time series. However, the effective implementation of some QRNN models is limited by the need of mid-circuit measurements. Those increase the requirements for quantum hardware, which in the current NISQ era does not allow reliable computations. Emulation arises as the main near-term alternative to explore the potential of QRNNs, but existing quantum emulators are not dedicated to circuits with multiple intermediate measurements. In this context, we design a specific emulation method that relies on density matrix formalism. The mathematical development is explicitly provided as a compact formulation by using tensor notation. It allows us to show how the present and past information from a time series is transmitted through the circuit, and how to reduce the computational cost in every time step of the emulated network. In addition, we derive the analytical gradient and the Hessian of the network outputs with respect to its trainable parameters, with an eye on gradient-based training and noisy outputs that would appear when using real quantum processors. We finally test the presented methods using a novel hardware-efficient ansatz and three diverse datasets that include univariate and multivariate time series. Our results show how QRNNs can make accurate predictions of future values by capturing non-trivial patterns of input series with different complexities.

## I Introduction

Processing and analysing multivariate time series, generally involve sophisticated algorithms that load high-dimensional datasets evolving in time, wherein temporal correlations are not trivial. In classical computation, Machine Learning is a consolidated approach for prediction and anomaly detection tasks [1; 2; 3; 4; 5; 6; 7]. Recurrent Neural Networks (RNNs) are a powerful tool for learning sequential data [8], and, over decades, several variations over the first model have improved their performance, like the well-known _Long Short-Term Memory_ (LSTM) or the _Gated Recurrent Unit_ (GRU) cells [1; 4; 9]. Transformers are a recent alternative [10] that seems to outperform the previous Neural Network models. However, there is not enough experience when using this algorithm for numerical multivariate time series.

The RNN models feature great results for sequential learning, but they are not problem-free. Some models cannot store information from first inputs in long time series. The LSTM cell arose to address this problem [1]. Moreover, because of the temporal correlations between different variables in multivariate time series, a high-dimensional space appears for computing data with non-linear patterns. Thus, it is necessary to build neural networks with more neurons, layers and parameters, which are computationally expensive and more challenging to train. Rapidly, we have to tackle a Deep Learning task, which requires computational resources and algorithms that ease the parameter optimisation during the neural network training, such as the backpropagation algorithm for estimating the gradients [11].

Quantum Machine Learning (QML) leverages the power of Quantum Computing to produce complex patterns that are probably not straightforwardly reproducible in a classical computer [12; 13]. However, in the current _Noisy Intermediate-Scale Quantum_ (NISQ) era of quantum computers, there is a need for algorithms requiring a low quantum circuit depth, since only a limited number of operations are meaningful. Consequently, a substantial part of the research in quantum algorithms is focusing on the Variational Quantum Algorithms (VQAs) [14]. A VQA lies on hybrid classical-quantum routines to classically minimise a cost function computed with a back-end Parameterised Quantum Circuit (PQC).

In this context, a VQA can be used as a Neural Network by feeding a quantum circuit with information from our dataset and then tuning the circuit parameters to minimise the cost function [15]. The first proposals of Quantum Recurrent Neural Networks (QRNNs) consist of a temporal loop that feeds a quantum state into a quantum circuit and applies the Schrodinger equation to evolve the state. Their applications were simulating stochastic processes [16] and stochastic filtering of signals with noise [17; 18]. In the middle of the NISQ era and with the explosion of Machine Learning, several proposals have arisen to enhance the power of the current Machine Learning models by adding Parameterised Quantum Circuits as part of the algorithm [19; 20; 21; 22; 23].

The actual QRNN model seeks to maximise the use and the power of Quantum Computing by a PQC that computes the whole sequence, and the only classical part
is data pre-processing, data post-processing and the optimisation of circuit parameters [24; 25; 26; 27; 28; 29]. The most common circuit architecture repeatedly encodes classical data into the quantum circuit, applies a unitary operator and measures a subset of qubits, involving multiple intermediate measurements before the end of the circuit. However, other proposals avoid intermediate measurements at the expense of increasing the number of qubits with the time-series length [27] or truncating the circuit [29].

The QRNN model is expected to provide advantages compared to classical neural networks, which extend to QML models in general. The encoding of classical data into quantum states allows us to compute a number of functions that increases exponentially with the number of qubits [15]. This feature enhances non-linearities that are required for learning complex patterns. In classical computing, the process requires exponential resources. Besides that, circuits with intermediate measurements, known as dynamic circuits, are not widespread in the QML literature, but they can lead to a new realm of algorithms on real quantum hardware [30], and they are being introduced in Quantum Computing platforms [31; 32].

The main aim of this article is to provide the mathematical tools to emulate this type of circuit, which is the core of QRNN algorithms. This is very interesting for two reasons: with the mathematical formulation, we can learn about the propagation of information through the quantum circuit, and the emulation makes the algorithm compatible with classical devices, before sending it to a quantum computer. Within this context, we test the theoretical procedure in several use cases, showing their potential for multivariate time series prediction.

The manuscript is structured as follows. In Section II we provide the method for emulating a quantum circuit with intermediate measurements, which has the property of returning values that depend on past data encoded in the circuit, like recurrent neural networks. As backpropagation is used in classical Machine Learning, in Section III we provide a method to analytically compute first- and second-order partial derivatives of the circuit outputs, by an expansion of the Parameter Shift Rule (PSR) [15; 33]. In Section IV we show the results after applying the former methods to an emulation of a QRNN for multivariate time series prediction, in order to demonstrate that the algorithm works by emulating the quantum circuit with the presented method. Finally, discussions are included in Section V.

## II Formulation of the QRNN states with density matrices

The QRNN model is based on the classical RNN model, as proposed by [25]. Following this approach, we explicitly provide a tensor representation of the internal states propagated through the QRNN and its outputs, based on the density matrix formalism. This representation allows the implementation of a classical emulator based on tensor operations and even simple matrix operations, which are very fast in current computational devices.

### The classical Recurrent Neural Network

Consider a set of data which is time-ordered,

\[\{\mathbf{x}_{(0)},\mathbf{x}_{(1)},\cdots,\mathbf{x}_{(t)},\cdots,\mathbf{x}_{(T)}\}, \tag{1}\]

that is the input multivariate time series, since each item is a vector containing the information of \(n_{v}\) variables. The output series is

\[\{\mathbf{y}_{(0)},\mathbf{y}_{(1)},\cdots,\mathbf{y}_{(t)},\cdots,\mathbf{y}_{(T)}\}, \tag{2}\]

which is the target in the Machine Learning task.

The output of a RNN at time \(t\) depends on the inputs from the previous time steps since it preserves past information with a form of memory [34]. The easiest way to imagine a network with memory is thinking in a box that continuously receives and returns data. To start, the box is supplied with an input, \(\mathbf{x}_{(0)}\), and it returns two objects: an output \(\overline{\mathbf{y}}_{(0)}\) which is read, and a _hidden state_\(\mathbf{h}_{(0)}\) which is re-introduced in the box next time. From the second time onwards, the box receives the previous hidden state \(\mathbf{h}_{(t-1)}\) and the input \(\mathbf{x}_{(t)}\). Then, it computes them and generates the output \(\overline{\mathbf{y}}_{(t)}\) and the hidden state \(\mathbf{h}_{(t)}\) to be re-introduced. The model is represented in Fig. 1, where we can see the recurrence, so that the behaviour is dynamic, in contrast to feedforward neural networks, for which the information is never transmitted backwards.

The information flux in a RNN through time splits into three different lines, represented in Fig. 1, that stand for _input data_\(\mathbf{x}_{(t)}\), _output data_\(\overline{\mathbf{y}}_{(t)}\) and _hidden state_\(\mathbf{h}_{(t)}\). Both \(\overline{\mathbf{y}}_{(t)}\) and \(\mathbf{h}_{(t)}\) are functions that depend on \(\mathbf{x}_{(t)}\) and \(\mathbf{h}_{(t-1)}\),

\[\left\{\begin{array}{ll}\overline{\mathbf{y}}_{(t)}&=\mathcal{Y}\left(\mathbf{x}_{(t )},\mathbf{h}_{(t-1)}\right)\\ \mathbf{h}_{(t)}&=\mathcal{S}\left(\mathbf{x}_{(t)},\mathbf{h}_{(t-1)}\right),\end{array}\right. \tag{3}\]

but at the same time, \(\mathbf{h}_{(t-1)}\) also depends on previous inputs and a hidden state. Functions \(\mathcal{Y}\) and \(\mathcal{S}\) are the part which is expensive to compute. They consist of matrix calculations that take into account the tunable parameters (weights), structure and connections between the layers of neurons that form the network.

### The Quantum Recurrent Neural Network

The QRNN leverages the power of PQCs implicitly performing matrix calculus. Like in the proposed classical model, wherein a network returns output data every time step, we measure the quantum circuit to obtain data
in our classical devices. Quantum Networks also contain a _hidden-state_ flux \(\mathbf{h}\). As its name suggests, a hidden state propagates internally and carries information not required to be known. That is why we can think of a quantum state as a hidden state. This quantum state is actually a mixed state that arises from the measurement of a subsystem A, and a subsystem B, which is not measured, provided that A and B are entangled.

Hence, we can now construct a general circuit that resembles the classical RNN described above, as it was proposed by [25]. Since a quantum circuit diagram represents the series of operations applied over several qubits (vertical direction) during the time (horizontal direction), the circuit in Fig. 2 can be interpreted as the unrolled representation through time of the QRNN. The circuit consists of two quantum registers: an _exchange register_ (E) with \(n_{E}\) qubits, and a _memory register_ (M) with \(n_{M}\) qubits. The former is used to exchange information between the quantum and the classical interface, by applying encoding and measurement operations, while the latter is never measured. The total number of qubits is \(n\). We then have a circuit _block_ for every time step: an estimation of an output from the input at time \(t\). The instructions are (see Fig. 2):

1. Initialise register E (M) to a desired pure state \(\ket{\Psi_{E}}\) (\(\ket{\Psi_{M}}\)). We will restrict the initialisation of all the qubits to \(\ket{0}\).
2. Reset of register E. All register E qubits start at \(\ket{0}\) every time step.
3. Apply a parameterised unitary (ansatz) \(U\left(\mathbf{x}_{(t)},\mathbf{\theta}\right)\) that evolves the state of each qubit and entangles some (or all) qubits from both registers E and M, correlating the information from both. \(\mathbf{\theta}\) are the variable parameters (weights) of the network.
4. Measure qubits from register E.
5. Repeat steps (2), (3), and (4) iteratively from \(t=1\) to \(t=T\).

The unitary \(U\left(\mathbf{x}_{(t)},\mathbf{\theta}\right)\) must encode the classical input \(\mathbf{x}_{(t)}\) each time step, apply parameterised gates depending on the set of weights \(\mathbf{\theta}\), and apply entanglement between registers E and M. The \(\mathbf{\theta}\) are always the same for the different repetitions of the unitary, following the classical approach (see Fig. 1).

Variations of this structure are available in the literature [28] when trying to create a circuit that sustains the coherence for a longer time. This is not the scope of this work, since we want to establish some bases for designing and emulating quantum RNNs. Solutions for NISQ-era limitations in these networks will require further research that may involve quantum hardware parameters, such as the readout features or the coherence itself [35].

### Memory and output of a QRNN

The way the vanilla RNN [36] (from which subsequent RNN algorithms were originated) operates data is the ideal sample of how data is explicitly processed in a RNN:

\[\begin{array}{ll}\mathbf{h}_{(t)}=&\sigma_{h}\left(W_{hx}\mathbf{x}_{(t)}+W_{hh}\bm {h}_{(t-1)}+\mathbf{b}_{h}\right),\\ \overline{\mathbf{y}}_{(t)}=&\sigma_{y}\left(W_{yh}\mathbf{h}_{(t)}+\mathbf{b}_{y}\right), \end{array} \tag{4}\]

being \(\{W_{hx},W_{hh},W_{yh}\}\) matrices with the weights, \(\{\mathbf{b}_{h},\mathbf{b}_{y}\}\) the biases, and \(\{\sigma_{h},\sigma_{y}\}\) the activation functions. Matrices (2-dimension tensors) are always represented by uppercase letters, while vectors (1-dimension tensors) are represented by bold lowercase letters.

A density matrix formalism is behind the classical emulation of the quantum circuit in Fig. 2. In this case, we can compute the exact probability distribution in every measurement of \(n_{E}\) qubits. At the same time, we include non-unitary operations, such as the measurement itself and the reset operation. An alternative is a state-vector emulation by sampling outputs, i.e., measurement probabilities are computed from multiple circuit executions, like in a real quantum device. However, exact probabilities are preferable when emulating circuits in VQAs, especially when the optimisation method is based on gradients.

Following the instructions to build the quantum circuits, we formulate the expressions to compute two ob

Figure 1: Classical Recurrent Neural Network representations. (a) Basic RNN scheme. (b) Basic RNN unrolled through time.
jects for every time step, as in Eq. (4): one is the reduced density matrix \(\rho_{M\,(t)}\), which transmits information from \(t\) to \(t+1\), while the other one is the expected value of some observable \(\left\langle O\right\rangle_{(t)}\) (the output before classical post-processing). Following the density matrix formalism, we get

\[\begin{array}{rl}\rho_{M\,(t)}=&\mathrm{Tr}_{E}\left[U\;\rho_{(t)}\,U^{ \dagger}\right],\\ \left\langle O\right\rangle_{(t)}=&\mathrm{Tr}\;\left[U\;\rho_{(t)}\,U^{ \dagger}\,O\otimes I^{\otimes n_{M}}\right],\end{array} \tag{5}\]

where \(\rho_{(t)}\) is the initial density matrix of the circuit at time \(t\) after the reset on register E and before applying the operator \(U=U(\mathbf{x}_{(t)},\mathbf{\theta})\). \(O\) is our observable, that will be, without loss of generality, diagonal. If the observable were not diagonal, we could get its spectral decomposition and apply the necessary transformations to the circuit(s) before measurement. The output (prediction) at time \(t\) is

\[\overline{y}_{(t)}=f(\left\langle O\right\rangle_{(t)}), \tag{6}\]

being \(f\) an arbitrary function. We are restricting to a single-variable output. However, the generalisation for multiple variables is straightforward by considering a set of multiple observables instead of a single one.

This section aims to derive a tensor representation of both \(\rho_{M\,(t)}\) and \(\left\langle O\right\rangle_{(t)}\) to (i) provide an explicit formula that can be implemented in a classical computer as matrix products, and (ii) show how the information is transmitted through the quantum circuit.

For the derivation of the formulas, a symbol and a group of \(r\) indices (lowercase letters) represent every mathematical object that is a \(r\)-rank tensor. The ordering of the indices follows the next criteria. A subindex in parenthesis refers to the current time step, and we can sometimes neglect it. The rest of the indices identify the coefficients for every projector or vector inside the Hilbert space. Indices above (below) correspond to the base vectors in the Hilbert space (dual Hilbert space). For quantum operators over the full circuit and density matrices representing the \(n\) qubits, indices are in pairs. The first index from the pair corresponds with register E, while the second one, with register M. We use the Einstein summation convention, i.e., terms are summed when upper and lower indices are repeated. Note that operations do not commute in general. See Appendix A for further explanation.

From now on, the dimension of the Hilbert spaces is represented by \(N_{E}=2^{n_{E}}\) and \(N_{M}=2^{n_{M}}\). The Hilbert space dimension corresponding to the \(n=n_{E}+n_{M}\) qubits is \(N=2^{n}\).

The reduced density matrix after measurement at a time \(t\) is

\[\begin{array}{rl}\left(\rho_{M}\right)^{m}_{(t)n}=&U^{im}_{(t)kq}\;\rho^{ kq}_{(t)lr}\;\left(U^{\dagger}\right)^{lr}_{(t)jn}\;\delta^{j}_{i}=\\ &\sum_{i=0}^{N_{E}-1}U^{im}_{(t)kq}\;\rho^{kq}_{(t)lr}\;\left(U^{ \dagger}\right)^{lr}_{(t)in}.\end{array} \tag{7}\]

By decomposing the operator into an encoding operator \(V(\mathbf{x}_{(t)})\) that acts only over register E and an operator \(W(\mathbf{\theta})\) that entangles both registers (see Fig. 3), we can separate the density matrix before applying \(W\) and then have that

\[\left(\rho_{M}\right)^{m}_{(t)n}=\sum_{i=0}^{N_{E}-1}W^{im}_{kq}\;\left(\left( \rho_{E}\right)^{k}_{(t)l}\left(\rho_{M}\right)^{q}_{(t-1)r}\right)\;\left(W^{ \dagger}\right)^{lr}_{in}, \tag{8}\]

where \(\left(\rho_{E}\right)^{k}_{(t)l}=V^{k}_{(t)0}\left(V^{\dagger}\right)^{0}_{(t )l}\). We have omitted \((t)\) in \(W\) since it does not depend on time. All the blocks include the same \(W\). We can see that \(\rho_{M\,(t)}\) depends on \(\rho_{M\,(t-1)}\) and \(\mathbf{x}_{(t)}\), as the hidden state in Eq. (4), unless the \(W\) operator is separable in the two registers (i.e. no entanglement between E and M).

Following Eq. (5), the expectation value of some diagonal observable \(O^{i}_{m}=d^{i}\delta^{i}_{m}\) (no sum

Figure 3: Decomposition of \(U\) operator into encoding (\(V\)) and evolution part (\(W\)). The latter is the same for all the circuit blocks since \(\mathbf{\theta}\) does not change during a circuit evaluation.

Figure 2: General form of the QRNN circuit. Arrows show the information flux.
symbol is above in both items) at time \(t\) is

\[\left\langle O\right\rangle_{(t)}=(\rho_{ij}^{\prime kl}O_{i}^{i}\beta_{ n}^{j})\delta_{kl}^{mn} \tag{9}\] \[=(\rho_{ij}^{\prime kl}d^{i}\delta_{m}^{i}\delta_{n}^{j})\delta_{kl} ^{mn}=(\rho_{ij}^{\prime kl}d^{i}\delta_{n}^{j})\delta_{kl}^{in}\]

where \(\rho^{\prime}\) is the density matrix after applying the \(U\) operator. By decomposing this density matrix, we have

\[\left\langle O\right\rangle_{(t)}=\sum_{i=0}^{N_{E}-1}d^{i}\sum_{n=0}^{N_{M}- 1}W_{kq}^{in}\left((\rho_{E})_{(t)l}^{k}(\rho_{M})_{(t-1)r}^{q}\right)\left(W^ {\dagger}\right)_{in}^{lr}, \tag{10}\]

which, again, proves the dependency of this observable with respect to both the inputs \(\mathbf{x}_{(t)}\) and the reduced density matrix from the previous step \(\rho_{M\;(t-1)}\), provided that \(W\) is an operator entangling E and M. In Eq. (4), it depends on the current hidden state, not the previous one. However, the recursion makes it dependent on \(\mathbf{h}_{(t-1)}\) too.

In most Quantum Machine Learning problems, the \(W\) operator is represented by a highly dense matrix (without null elements) and difficult to decompose since in general it should be formed by several entanglement layers. Then, during emulation, it is plausible to build this operator, keep it in memory and then we can split it into parts and operate with the parts as matrix products, which are very fast operations in classical computers. That splitting can be done by creating \(W^{i}\) matrices, which are \(N_{M}\times N\) dense matrices. With this derivation, we do not need to compute the complete \(N\times N\) density matrix before computing \(\rho_{M\;(t)}\) and \(\left\langle O\right\rangle_{(t)}\). Moreover, these operations are performed only once for both Eqs. (8) and (10), because the sum over \(n\) in Eq. (10) extracts the items with \(m=n\) from Eq. (8). In a matrix representation, this is the trace of each of the \(N_{E}\) matrices that summed together give rise to \(\rho_{M\;(t)}\).

## III Analytical derivatives

In Machine Learning, most parameter optimisation algorithms are based on gradients, which require some method to be computed. One is numerical differentiation, which is very inaccurate, and another one is symbolic differentiation, which is more computationally expensive [37]. Automatic Differentiation, with _backpropagation_ as the most-known algorithm, seems to have many advantages.

In QML, as current circuits are prone to errors and the outputs are obtained after running the circuit multiple times (shots), systematic and stochastic errors seriously affect the estimation of the cost functions. Errors are amplified when computing numerical gradients, which makes it much more difficult to compute them unless we run circuits with a great number of shots and apply sophisticated error mitigation techniques. Automatic Differentiation as is done in classical computing is not possible, because we need to store intermediate results through the network and intermediate quantum states cannot be accessed in real quantum devices without affecting them. To access them, we need to measure them, and they collapse.

Despite these restrictions, the Parameter Shift Rule (PSR) [15; 33] provides a method to compute analytical gradients in quantum hardware. Hereafter, we consider circuits with parameters encoded into rotation gates \(R_{i}(\alpha)=e^{-i\alpha\sigma_{i}}\) generated by a Pauli matrix \(\sigma_{i}\), because they are a set of fundamental gates used for making PQCs in gate-based quantum devices.

### First-order partial derivatives

We define the shift expectation value at time \(t\), \(\left\langle O\right\rangle_{(t)}^{\chi}\Big{|}_{ri}\), as the expectation value at time \(t\) after shifting the parameter \(\theta_{i}\) a value \(\chi\) at block \(t_{r}<t\) in the quantum circuit. We recall that if \(t_{r}>t\), the shift does not affect the output at time \(t\).

For the set of gates that we consider, the shifts are \(\pm\frac{\pi}{2}\)[33]. Then, in the remaining, we write only the sign.

The partial derivative of the observable at time \(t\), \(\left\langle O\right\rangle_{(t)}\), is then

\[\partial_{i}\left\langle O\right\rangle_{(t)}\equiv\frac{\partial\left\langle O \right\rangle_{(t)}}{\partial\theta_{i}}=\sum_{r=0}^{t}\frac{1}{2}\left.\left( \left\langle O\right\rangle_{(t)}^{+}-\left\langle O\right\rangle_{(t)}^{-} \right)\right|_{ri}, \tag{11}\]

derived in Appendix B.

In classical RNNs, the chain rule propagates backwards in time, giving rise to _backpropagation through time_ (BPTT) [38]. This method is necessary for computing gradients in RNNs, as backpropagation is used in deep learning models, but it has some caveats. Firstly, it requires computing plenty of terms because we need to repeat the network the same number of time steps we have. When contributions from the beginning of the series are negligible, a possible solution is truncation methods [39]. Secondly, the propagation backwards in time requires several weights-matrix products. Instability manifests as _vanishing_ (_exploding_) gradients if the eigenvalues of the weights matrices are lower (higher) than \(1\)[40]. Truncation can partially address this problem [41; 39].

From equation (11), computing the partial derivative of the observable at time \(t\) requires \(2t\) function evaluations with this method. Therefore, computing the gradient requires \(2tN_{\theta}\) function evaluations, being \(N_{\theta}\) the number of parameters of our ansatz. The loss function for training depends on the \(T\) time steps. We can evaluate the outputs in a single circuit evaluation, i.e. running the circuit a given number of times (shots) but fixing all the parameters. Thus, in general, the number of function evaluations needed to compute the exact gradient of the loss function is \(2TN_{\theta}\). This is worse than the case of circuits without mid-circuit measurements. Previous proposals implement techniques to reduce the number of
circuit evaluations [42; 43; 33], but not for this type of intermediate-measurement-based circuits.

Instability problems cannot appear in Quantum Neural Networks in the form of _exploding gradients_ because Quantum Circuits naturally implement unitary operations [44] before measuring. Unitary operations are a restriction of some classical ML models, precisely to avoid exploding gradients, and there are ways to implement them without losing too much expressivity [45]. In the case of a QRNN, it is straightforward to see that exploding gradients cannot appear, by looking at equation (11) since the subtraction of two observables cannot blow up, and the partial derivative only involves the summation of these terms.

The _vanishing gradient_ problem does appear in VQAs, where parameterised circuits can suffer from barren plateaus in the optimisation landscape [46; 47]. In the case of RNNs, vanishing gradients cause loss of memory: the network stops depending on past inputs at some moment during the training. In the QRNN model, dependencies of past inputs gradually vanish too, because, at time \(t\), the network starts in a mixed state, \(\rho_{E\,(t)}\otimes\rho_{M\,(t-1)}\), evolved by a unitary operator. The terms from \(\rho_{M\,(t-1)}\) are attenuated after these operations and many circuit blocks.

### Second-order partial derivatives

Some optimisers require the computation of the Hessian matrix, apart from the Jacobian (gradient). Some authors use the PSR to compute the Hessian [48; 49] in PQCs. Here, we provide the method for computing second-order partial derivatives in the QRNN, by applying the PSR.

We define the double-shift expectation value at the time \(t\) as \(\left\langle O\right\rangle_{(t)}^{\chi\lambda}\right|_{sj}^{ri}\), which is the expectation value at time \(t\) after shifting the parameter \(\theta_{i}\) in block \(t_{r}\) a value \(\chi\) and the parameter \(\theta_{j}\) a value \(\lambda\) in the block \(t_{s}\) of the quantum circuit.

The second-order partial derivatives are

\[\partial_{i}\partial_{j}\left\langle O\right\rangle_{(t)}\equiv \frac{\partial^{2}\left\langle O\right\rangle_{(t)}}{\partial\theta_{i}\partial \theta_{j}} \tag{12}\] \[-\frac{1}{4}\sum_{r}^{t}\sum_{s}^{t}\left.\left(\left\langle O \right\rangle_{(t)}^{+-}+\left\langle O\right\rangle_{(t)}^{-+}\right)\right| _{sj}^{ri},\]

for \(i\neq j\).

Moreover,

\[\partial_{i}^{2}\left\langle O\right\rangle_{(t)}\equiv\frac{ \partial^{2}\left\langle O\right\rangle_{(t)}}{\partial\theta_{i}^{2}} \tag{13}\] \[-\frac{1}{2}\sum_{r}^{t}\sum_{s}^{r-1}\left.\left(\left\langle O \right\rangle_{(t)}^{+-}+\left\langle O\right\rangle_{(t)}^{-+}\right)\right| _{si}^{ri}-\frac{1}{2}t\left\langle O\right\rangle_{(t)}\]

for \(i=j\). These are derived in Appendix B.

Considering these expressions, the symmetry \(\partial_{i}\partial_{j}=\partial_{j}\partial_{i}\) and that we can recycle the intermediate calculations for the final estimation of the Hessian of the loss function, the total number of function evaluations is \(2T^{2}N_{\theta}^{2}+1\) (probe in Appendix B). The symmetries allow us to reduce the computational cost. However, the scaling with time steps and number of parameters is quadratic. Some approximations would be needed to improve this scaling.

Both gradient and Hessian calculations are easy to parallelise. Then, it would be plausible to simultaneously use multiple Quantum Processing Units (QPUs) to perform several circuit evaluations. The problem is the high noise level of the current quantum hardware and different noise models for different QPUs that could complicate the optimisation process. Further investigation is required to see the scope of these architectures.

## IV Results

We have implemented an algorithm that emulates a Quantum Recurrent Neural Network, by using the definitions and methods described in the previous sections. The model is based on a quantum circuit that uses a hardware-efficient ansatz. This ansatz uses layers of rotation gates for encoding classical inputs into the quantum circuit, and alternates layers of single-qubit-rotation gates parameterised by the set \(\mathbf{\theta}\) and CZ gates.

The explicit circuit ansatz \(U(\mathbf{x}_{(t)},\mathbf{\theta})\), used to test the performance of the general model, is represented in Fig. 4. The encoding part \(V(\mathbf{x}_{(t)})\) acts only over register E qubits and repeats the encoding for each input value \(R+1\) times in order to have a better expressivity [50; 51; 52]. The evolution and entanglement operator \(W(\mathbf{\theta})\) does not vary during the circuit since it does depend always on the same parameters. We try to maximise its expressibility by repeating the entangling layers several times with different parameters.

In order to test the ideal QRNN emulation, and its learning capabilities, we use three datasets: (a) a dimmed triangular signal, (b) a non-linear damping signal with a sinusoidal perturbation, and a set (c) consisting of two non-linear damping signals as input and a linear combination of them as output. These three examples are meaningful for two reasons: the signals are nonlinear,
a feature that hinders the training, and are different enough to test several quantum circuit configurations. Moreover, the third one demonstrates the applicability of the model to, at least, two-variable series. Details about data generation are included in Appendix C.

The network is trained by optimising the set of parameters \(\mathbf{\theta}\), using the L-BFGS-B algorithm [53]. We divide the series into windows of \(T=20\) points and the task is predicting 5 points for the output variable. The windows are divided into three sets: 20 % for testing, 20 % of the remaining for validation, and the rest for training. The distribution of validation samples is the same for the three datasets, but randomly generated. The loss function is the Root Mean Square Error (RMSE) between the output \(\overline{y}\) and the target series \(y\) of all the training windows. The prediction for every time step is

\[\overline{y}_{(t)}=\left\langle Z^{\otimes n_{E}}\right\rangle_{(t)}+b, \tag{14}\]

where \(Z\) is the \(\sigma_{z}\) Pauli matrix and \(b\) is a bias, which is an extra trainable parameter. These values are computed from an exact expectation value, without statistical variations typically arising when sampling the outputs in a real quantum computer.

The network hyperparameters used for each case are indicated in Table 1. They were not optimised during this analysis, since the scope of this section is only to test the algorithm. Further works would assess different ansatz architectures, hyperparameters and optimisation techniques.

We have executed 8 different parameter optimisations for each dataset, randomly initialising the rotation gates parameters in the interval \([0,1)\). The set of parameters after optimisation strongly depends on the initialisation. That has consequences in the training, leading to differences in the prediction accuracy. Nonetheless, results never show a relative RMSE greater than 10 % of the series range (from -0.75 to 0.75). The set of parameters leading to the lower RMSE in the validation set are selected as the solution among the 8 different optimisation results. The corresponding datasets and their predictions are shown in Fig. 5. Table 2 contains the resulting RMSE of the estimated points with respect to the training, validation and test targets. The fourth RMSE corresponds to the prediction of the whole test series, not only the last 5 points of each window, by shifting the input windows 5 by 5 points.

The L-BFGS-B optimiser uses the gradient of the loss function to find its minimum. This algorithm does not compute the complete Hessian, but it makes an approximation. Then, the analytical Hessian is not needed. Furthermore, the analytical gradient can be substituted by a numerical one if the precision of the loss function is high enough. We use the default 2-point finite difference estimation with an absolute step size \(\epsilon=1\times 10^{-8}\) implemented in the optimiser, as well as the analytical form of the gradient. For the latter, we evaluate the partial derivatives of the circuit outputs by the formulas provided in Section III, and then we apply the chain rule to calculate the gradient of the loss function.

We have compared between analytical and numerical optimisation, as shown in Fig. 6 and Table 2. The op

\begin{table}
\begin{tabular}{c|c|c|c|c|c|} Case & \(n_{E}\) & \(n_{M}\) & \(L\) & \(R\) & \(N_{\mathbf{\theta}}\) \\ \hline (a) & 1 & 2 & 2 & 3 & 31 \\ (b) & 2 & 3 & 2 & 1 & 43 \\ (c) & 2 & 3 & 5 & 3 & 100 \\ \end{tabular}
\end{table}
Table 1: Quantum circuit configuration for each case analysed. In case (b), we re-upload input data in two qubits.

Figure 4: QRNN ansatz, \(U(\mathbf{x}_{(t)},\mathbf{\theta})\), consisting of two parts. The first one is the data encoding, and gates inside the orange box are repeated with different parameters, that are a subset of trainable parameters, \(\mathbf{\alpha}_{i}^{r}\in\{\mathbf{\theta}\}\). We use use one qubit per input variable. The second one is the evolution and entanglement part, where the blue box is repeated \(L\) times (layers). Each layer is a column of \(U_{3}\) rotations parameterised by a threesome of parameters, \(\mathbf{\beta}_{i}^{l}\in\{\mathbf{\theta}\}\), and CZ gates entangling every qubit from E with every qubit from M. A final column of \(U_{3}\) gates is applied over register E before measurement.
timisation converges to close RMSE values for both the training and validation loss curves, despite following a different trajectory in some zones. Although analytical calculations are commonly used in neural networks to improve accuracy [37], we cannot ensure that they are better in the QML tasks provided here, when emulating the circuit without any noise. Apart from that, the plots show a quick convergence in case (a), but slower in cases (b) and (c). Validation curves (and final loss values) are above the training curves in cases (a) and (c). However, they remain relatively stable and close to the training ones. Case (b) manifests validation RMSE values below the training ones, due to the distribution of validation samples. A different election would return distinct validation values. In fact, a higher number of samples should reduce the dependence on the election of validation samples. In contrast, case (c) manifests a training RMSE that slightly moves downwards from the validation one, due to a little overfitting that is not meaningful in this context. Those issues should be addressed in a complete machine learning project. Regarding they are beyond the aim of this work, we will leave them for future research.

In ideal emulation without sampling (exact), the emulated quantum circuit returns an expectation value with a precision given by the classical machine variables (e.g. float-64). As, in general, an expectation value cannot be extracted from a quantum circuit in a single shot, its evaluation by multiple repetitions (shots) of the circuit leads to stochastic variations due to this sampling process. That makes numerical gradients very inaccurate and we must use analytical gradients. The PSR allows us to tackle this problem in real quantum circuits, but

\begin{table}
\begin{tabular}{c c|c|c|c|c|c|c|c} \multirow{2}{*}{**Case Gradient**} & \multirow{2}{*}{\(n_{\textbf{it}}\)} & \multirow{2}{*}{\(n_{\textbf{lev}}\)} & \multirow{2}{*}{\(n_{\textbf{lev}}\)} & \multirow{2}{*}{**RMSE tr.**} & \multirow{2}{*}{**RMSE val.**} & \multirow{2}{*}{**RMSE tes.**} & \multirow{2}{*}{**RMSE fte.**} \\ \cline{4-9}  & & & & & & & & \\ \cline{4-9}  & \multirow{2}{*}{analytical} & 600 & 701 & 701 & 0.006 & 0.011 & 0.004 & 0.005 \\ \cline{4-9}  & & 553 & 20896 & 653 & 0.005 & 0.008 & 0.003 & 0.004 \\ \cline{4-9} \multirow{2}{*}{(b)} & \multirow{2}{*}{analytical} & 579 & 664 & 664 & 0.091 & 0.038 & 0.118 & 0.082 \\ \cline{4-9}  & & 581 & 28556 & 649 & 0.090 & 0.041 & 0.125 & 0.082 \\ \cline{4-9} \multirow{2}{*}{(c)} & \multirow{2}{*}{analytical} & 1000 & 1082 & 1082 & 0.035 & 0.048 & 0.050 & 0.044 \\ \cline{4-9}  & & 1000 & 110393 & 1082 & 0.036 & 0.048 & 0.051 & 0.045 \\ \hline \end{tabular}
\end{table}
Table 2: Optimisation solutions for each case with analytical and numerical gradients: number of iterations, number of function evaluations, number of Jacobian evaluations, and RMSE in training, validation and (complete) test sequences.

Figure 5: Results of the learning task with each one of the three datasets. Here, numerical gradients were used. Each series was divided into windows of 20 points, each one being a sample. The neural network must predict the value of the last 5 outputs for each window, by reducing the RMSE with respect to the last 5 points of the target. Dashed lines are the inputs \(\mathbf{x}_{(t)}\), from which we make the prediction. Solid lines represent the targets. Points are the predictions. Windows are represented in orange (red), blue (dark blue) or green (dark green, pink) depending on whether the window is used for training (tra), validation (val) or testing (tes), respectively. In the test region, we add a window to make the prediction for the full test sequence (fte).
ideal exact emulation using numerical gradients requires less computational resources for the cases studied. The density matrix emulation enables exact probabilities calculation, but also simulating stochastic noise to investigate its effect on the algorithm's performance. In the latter case, we would need analytical gradients.

The final loss values and the plots in Fig. 5 show the ability of the network to model non-trivial patterns in univariate and multivariate input series. The univariate output series estimations approximate to the output targets under a reasonable error because the RMSE is, in the worst case, an order of magnitude lower than the range of output values (from -0.75 to 0.75). The lowest error corresponds with the univariate triangular signal, because it is the simplest case. Note that, in this case, the test RMSE is lower because the series decreases in amplitude with time. Meanwhile, the errors for the damped oscillator signals increase to upper values, but the network can reliably extrapolate the time series to the test region, for both the 1-variable and 2-variable cases. Apart from that, the iterations that the optimiser needs during the training are comparable. Furthermore, the more complex, the more parameters are needed, so that qubits and layers are added with increasing complexity from the case (a) to (c). In short, our results confirm the actual applicability of QRNN models to univariate and multivariate time series prediction and their emulation through the density matrix formalism.

## V Discussion

We have used the density matrix formalism to derive the formulation of the hidden states and the outputs in a Quantum Recurrent Neural Network with intermediate measurements, and its first- and second-order partial derivatives with respect to its trainable parameters. With the density matrix formulation, it is possible to directly perform an ideal emulation of the parameterised quantum circuit that makes up the network and we showed the results for a machine learning task with three different datasets. The analytical gradients will permit to train networks with noisy outputs coming from quantum circuits.

The QRNN was presented in previous works [24, 25, 26, 27, 28, 29]. Some of them delve into the inner structure of such quantum circuits, but there was not an explicit formulation to see how quantum states are operated through the quantum circuit. With our work, we contribute to a better understanding of how information is processed inside the network, and we compare it with the classical RNN. Moreover, the given formulas allow a direct implementation of a code for emulation.

The formulation provides various remarkable advantages related to other possible methods. First of all, it avoids computing the complete density matrix of the \(n\)-qubits quantum state, whose number of elements scales as \(\mathcal{O}(4^{n})\), by splitting the ansatz quantum operator. Moreover, the network outputs are directly computed from the terms that build up the reduced density matrix. Secondly, density matrix formulation allows the computation of exact expectation values in circuits with intermediate measurements, which is not feasible with state-vector emulation for real cases of interest.

Consequently, the gradients can be computed with high precision in a classical computer. We have seen that parameter optimisations with numerical and analytical gradients converge to close values. The great difference between both methods is the computational cost. Analytical forms require \(2T\) circuit evaluations for a gradient component, being \(T\) the number of time steps (circuit blocks inside a single circuit). Meanwhile, a 2-point finite difference approximation simultaneously shifts the parameter in all circuit blocks, needing only 2 circuit evaluations per gradient component. Nonetheless, as the precision of the variables is lost, analytical forms become necessary. Then, we expect that they will be unavoidable in emulations with noise and executions in real quantum computers. Our formulation lets us use numerical gra

Figure 6: Curves of optimisation for each of the three datasets, showing how the value of the RMSE varies with the iterations, for training and validation sets, and comparing results with both analytical and numerical gradients. The optimisations are those with lower final validation RMSE.
dients due to its high precision, at least for the datasets that were studied. Future work may explore the effect of quantum circuit noise. Similar conclusions apply to the Hessian calculation, which requires a \(T^{2}\) factor of circuit evaluations to be computed analytically, but it is usually approximated due to its high computational cost in most cases.

In this context, the idea of distributed quantum computing, where several Quantum Processing Units work simultaneously, gathers strength. Despite the number of circuit evaluations needed, the QRNN model shows a significant feature: it does not require a lot of qubits. For the hardware-efficient ansatz and the datasets we studied, an ideal quantum processor of 5 qubits would be enough to run the circuits.

The results are promising for potential use in real cases. To this end, the model hyperparameters optimisation to enhance the power of the neural network, the adaptability to the quantum hardware limitations, research on circuit structures that capture correlations behind complex datasets and the adaptability of different optimisation techniques are starting points towards a generalised application of Quantum Recurrent Neural Networks to multivariate time series prediction.

###### Acknowledgements.

 We thank the CESGA Quantum Computing group members for their feedback and the stimulating intellectual environment they provide. We thank Constantino Rodriguez Ramos, especially for theoretical insights. This work was supported by Axencia Galega de Innovacion through the Grant Agreement "Despegamento dunha infraestructura baseada en tecnodoxias cuanticas da informacion que permitia physarum a I+D+I en Galicia" within the program FEDER Galicia 2014-2020. A. Gomez was supported by MICIN through the European Union NextGenerationEU recovery plan (PRTR-C17.11), and by the Galician Regional Government through the "Planes Complementarios de I+D+I con las Comunidades Autonomas" in Quantum Communication. Simulations on this work were performed using the Finisterrae III Supercomputer, funded by the project CESGA-01 FINISTERRAE III.

## Author Contributions

J. D. V., M. M. J. and A. G. conceived the problem. J. D.V. developed the mathematical model and programmed the emulator. D.F. reviewed the mathematical model. All the authors contributed to analyse the data and review the manuscript.
