1 Motivation and introduction

Learning quantum dynamics presents fundamental challenges in quantum physics, to which numerous machine learning techniques have been applied (Huang et al. 2023; Caro et al. 2023; Flurin et al. 2020; Carleo and Troyer 2017; Jerbi et al. 2023; Mohseni et al. 2022; Granade et al. 2012; Rodríguez et al. 2022, and the references therein). Significant efforts have been made to compactly represent many-body quantum states and efficiently parameterize their dynamics using classical learning algorithms with classical data (Carleo and Troyer 2017; Mohseni et al. 2022; Schmitt and Heyl 2020; Hartmann and Carleo 2019). On the other hand, quantum algorithms that work on a large batch of classical data typically require extreme assumptions about data loading and readout, which remain points of contention, see, for example, Aaronson (2015); Tang (2019); Matteo et al. (2020).

The two challenges spur the proposals of quantum algorithms that can learn directly from quantum data (Schatzki et al. 2021; Huang et al. 2022). Our work adds to the literature by providing a novel quantum algorithm for learning from quantum data inspired by next generation reservoir computing (NG-RC) algorithm. Importantly, this model-free paradigm of reservoir computing requires time-series of quantum data, without assuming a quantum dynamical ansatz a priori. Such approach circumvents inductive biases that could arise in learning complex quantum dynamics. We commence with a concise overview of the standard RC and the NG-RC, emphasizing the innovative utilization of NG-RC in the accurate prediction of many-body quantum dynamics. Subsequently, we address the limitations inherent to classical NG-RC and introduce a quantum algorithm designed to overcome the challenges presented by its classical counterpart.

1.1 Reservoir computing (RC)

RC is a computational paradigm in machine learning that harnesses a recurrent neural network (RNN) to learn time-series data, such as the states of a dynamical system. Even when the dynamics is complicated or chaotic, well-tuned RC can accurately forecast the future states of such dynamics up to many Lyapunov times (Pathak et al. 2017, 2018). Importantly, RC bypasses the training of an RNN by utilizing a fixed, randomly initialized RNN, called a reservoir, consisting of a large number L of hidden neurons.

Suppose that an input data \(\varvec{s}_k\in \mathbb {R}^{D}\) is fed into the reservoir, where k specifies the time step of some dynamical process. The data is represented as a feature vector \(\varvec{x}_k = f(\varvec{s}_k,\varvec{x}_{k-1}) \in \mathbb {R}^L\), where f is typically a highly non-linear function that represents the dynamics of the reservoir. The feature vector is then linearly transformed in the final, trainable output layer into a prediction vector \(\varvec{y}'_k = W\varvec{x}_k\). In particular, the input \(\varvec{s}_k\) at the kth time step and the output at the previous time step \(\varvec{y}'_{k-1}\) may be fed together as inputs into the reservoir to train the feature vector \(\varvec{x}_k\) for the next time step.

To prevent overfitting in a supervised learning via RC, the weight matrix W is obtained by optimizing \(\varvec{y}'_k\) in the least-square sense with respect to the desired target \(\varvec{y}_k\), which can be done via the Tikhonov regularization

$$\begin{aligned} W = YX^T(XX^T+\lambda I)^{-1}. \end{aligned}$$
(1)

Here, \(Y=(\varvec{y}_0,\ldots ,\varvec{y}_{T-1})\in \mathbb {R}^{D\times T}\) is the target matrix, \(X = (\varvec{x}_0,\ldots ,\varvec{x}_{T-1})\in \mathbb {R}^{L\times T}\) is the feature matrix, and \(\lambda \ge 0\) is the regularization parameter. Here, T stands for the number of time steps of dynamics used to train the reservoir.

1.2 NG-RC

Despite the fast training protocol offered by RC, the random initialization of the reservoir presents its own problem: there is an overwhelmingly large number of hyperparameters to be optimized, and there is no consensus on how to pick an optimal reservoir. NG-RC is an alternative approach that takes advantage of the discovery that RC can be equivalently performed using a linear reservoir and a nonlinear trainable output layer (Gonon and Ortega 2020; Hart et al. 2021). The latter is in turn equivalent to a nonlinear vector autoregression (NVAR) machine (Bollt 2021). An NVAR machine predicts future observations of a time series using past observations. In particular, the underlying state space of the dynamics to be learned can be reconstructed using linear and nonlinear functions of past observations. Inspired by this correspondence, NG-RC (Gauthier et al. 2021) proposes taking m time-delay data:

$$\begin{aligned} \varvec{o}_k=\varvec{s}_{k}\oplus \varvec{s}_{k-\Delta }\oplus \varvec{s}_{k-2\Delta }\oplus \ldots \oplus \varvec{s}_{k-(m-1)\Delta } \end{aligned}$$
(2)

with step size \(\Delta \) as an input, and, forgoing the need for a reservoir, directly constructing a feature vector,

$$\begin{aligned} \varvec{x}_k=\varvec{o}_k\oplus (\varvec{o}_k)^{\otimes p}, \end{aligned}$$
(3)

whose nonlinearity arises from a degree-p monomial of the m-delay data, \(\Delta \), m, and p being the NG-RC hyperparameters. The feature vector is then optimized via a linear least-square regularization to predict the target dynamics.

NG-RC with \(\Delta =1\), \(m=2\), and \(p=2\) has been used to efficiently predict the dynamics of the Lorenz attractor using only small data sets (Gauthier et al. 2021). In fact, with this set of hyperparameters, which possess minimal memory of the past and a minimal degree of non-linearity, the NG-RC can still perform very well, as we demonstrate later in this paper, in forecasting the dynamics of both integrable and chaotic Ising models. The feature vector in this case is a \((4D^2+2D)\)-dimensional vector of the form

$$\begin{aligned} \varvec{x}_k = \varvec{s}_{k}\oplus \varvec{s}_{k-1} \oplus [(\varvec{s}_{k}\oplus \varvec{s}_{k-1}) \otimes (\varvec{s}_{k}\oplus \varvec{s}_{k-1})]. \end{aligned}$$
(4)

While NG-RC has been demonstrated to effectively forecast complex classical dynamics such as spatiotemporal chaos (Barbosa and Gauthier 2022) and exhibits robust predictive capabilities despite noisy input dynamics (Liu et al. 2023), its potential application to forecasting the dynamics of quantum observables and more generally of full many-body quantum states, which encompass both real and imaginary components, remains unexplored.

1.3 Learning quantum dynamics with NG-RC

Here, we demonstrate the capability of NG-RC to predict a unitary quantum dynamics using the strategy we refer to as skipping ahead. Instead of training the output of the feature vector Eq. 4 to match the target state at the next time step, \(\varvec{y}_k = \varvec{s}_{k+1}\), we could train it to match the target state at the next \(\tau \) steps, skipping into the future. That is, given \(\{\varvec{x}_k\}_{k=0}^{T-1}\) according to Eq. 4, the next generation reservoir is trained using \(\{\varvec{y}_k = \varvec{s}_{k+\tau }\}_{k=0}^{T-1}\). We found that our approach of skipping ahead yields a remarkably accurate prediction of the full quantum state into the far future, retaining the fidelity measure close to 1 in the standard example of integrable one-dimensional many-body dynamics that we now discuss.

The one-dimensional transverse-field Ising model is governed by the Hamiltonian

$$\begin{aligned} \mathcal {H}=-J\sum _{i=1}^d Z_{i}Z_{i+1} + h\sum _{i=1}^d X_{i}, \end{aligned}$$

with a periodic boundary condition, where \(X_i\) (resp. \(Z_i\)) is the Pauli-X (resp. Z) operator acting on the ith qubit. The system can be exactly solved by mapping to the free-fermionic model via the Jordan-Wigner transformation. However, the computation of certain expectation values such as the two-point correlation function \(\langle {X_iX_j}\rangle \), let alone the computation of the time-evolved state, on a classical computer becomes infeasible as the number of qubits becomes large (Bravyi and Kitaev 2002; Beenakker et al. 2004). (Note that the model can be simulated efficiently on a quantum computer by exact diagonalization (Cervera-Lierta 2018).)

We benchmark the prediction for the time evolution of a four-qubit transverse-field Ising model in the disordered phase with \(J=0.5\) and \(h=5\). The initial quantum state \(|{0000}\rangle \) followed by a long burn-in period initializes the inputs for the predictions of unseen \(\widetilde{T}=4\times 10^4\) time steps in the future. (After the burn-in period, the input states are far into the future of those from the training set.) Fig. 2 (left column) shows the results employing the skipping-ahead strategy with the time skip of \(\tau =10^6\) steps into the future, whereas Fig. 2 (right column) shows the results from the conventional iterative NG-RC approach to predict successive states in the future. Notably, the skipping-ahead strategy significantly outperforms the conventional approach as reflected by the almost perfect fidelity between the predicted states and the target states, as well as the almost perfect prediction of the observables. The iterative approach’s performances drop as the number of iteration increases, as expected.

Fig. 1
figure 1

(Top) The schematic of the QNG-RC algorithm. The initial part of the algorithm encodes the history of quantum dynamics collected from a time series into nonlinear feature vectors \(|x_k\rangle \), which is then processed by a regularized linear layer to predict future quantum states. (Middle) In the training phase of QNG-RC, the optimally trained weight matrix in Eq. 1 is encoded in a unitary operator and processed by the quantum singular value transform (QSVT), the total complexity of which is summarized in Theorem 2. (Bottom) In the prediction phase, the feature vector can be constructed via the same quantum circuits as in the training phase. Consequently, the predicted states will be revealed in the blue register in the last quantum circuit after applying the optimal weight matrix from the training phase

Full size image

In addition, we also benchmark our skipping-ahead approach with a non-integrable, chaotic many-body dynamics in a tilted-field Ising model with five qubits. The fidelity of the predicted states remains above 0.99, showing an excellent agreement with the full many-body states into the future. The details and results are provided in Section S5.

1.4 Challenges in using NG-RC for quantum dynamics predictions

The preceding results provide a piece of evidence that NG-RC may offer an alternative route to a high-precision prediction of quantum dynamics. However, the primary computational bottleneck arises from matrix inversion in Eq. 1 (or its generalization for complex-valued data with \(X^T\) replaced by \(X^{\dagger }\)). The Tikhonov-regularized least squares require \(O(MN^2)\) matrix operations, where N and M are, respectively, the larger and the smaller dimension of the matrix. This complexity is computationally prohibitive when the data matrix is collected from many-body quantum states, since the dimension of each column vector scales exponentially with the system size. In the following section, we present an end-to-end, quantum next generation reservoir computing (QNG-RC) algorithm that does not suffer from the exponential complexity of the classical counterpart.

1.5 Related works

A significant portion of work on quantum variants of classical reservoir computing, referred to as quantum reservoir computing (QRC), utilizes quantum systems as reservoirs (substrates) (Fujii and Nakajima 2017; Martínez-Peña et al. 2024; Nokkala et al. 2021; Mujal et al. 2021; Fujii and Nakajima 2021; Kutvonen et al. 2020; Martínez-Peña et al. 2021; Mujal et al. 2023; Martínez-Peña and Ortega 2023; Kubota et al. 2023; Fry et al. 2023). However, the outputs of these QRCs remain classical in nature, namely the expectation values of some observables, which are subsequently used to predict classical dynamical systems, not the full quantum state dynamics as in this work. Other QRC proposals focus on the task to extract partial information about quantum states such as their entanglement or purity (Ghosh et al. 2019) or perform a complete tomography (Ghosh et al. 2020). Kawai and Nakagawa (2020) propose learning excited states of a Hamiltonian by measuring expectation values in the ground state, which are then subsequently optimized using classical machine learning techniques. In all three examples, the main goal is to use QRC to effectively transform final single-qubit measurements into entangled ones that contain the desired information about the quantum state of interest, thus reducing the number or the complexity of measurement settings.

Nevertheless, there exist notable works that share several characteristics with ours. Carleo and Troyer (2017) introduce a method for learning quantum states and dynamics using neural network states, circumventing inductive biases and the exponential size of classical memory needed to process quantum data. Since reservoirs can be considered as special cases of recurrent neural networks, conventional QRC cannot surpass the predictions of the best quantum neural networks. However, the training of neural networks itself is a complex task in its own right, and both RC and NG-RC offer a means to navigate this challenge. Moreover, Ref. Carleo and Troyer (2017) learn a quantum dynamics by classically processing the data obtained by sampling from the trained network, a point that diverges from the approach presented in our work. Another related contribution, highlighted in Ghosh et al. (2021) employs QRC to coherently simulate the target state of a quantum dynamics, whereas our approach shifts the methodology to that of a next generation RC.

While our skipping-ahead method offers a fast route to generate quantum states that would otherwise emerge in the far future through the natural progression of quantum dynamics, the technique should not be confused with the fast-forwarding of quantum evolution, since our method requires training samples from the equally distant past. Moreover, our prediction operates in a model-free framework, devoid of Trotterization techniques. This distinction ensures that our proposed skipping-ahead approach is fundamentally different from fast forwarding. Consequently, the constraints outlined in the no-fast-forwarding literature (Cîrstoiu et al. 2020; Berry et al. 2007; Atia and Aharonov 2017; Haah et al. 2018; Gu et al. 2021) should not apply to our prediction method. Additionally, unlike many QRC proposals that utilize near-term quantum devices as quantum reservoirs, our circuit-based QNG-RC is designed for implementation on fault-tolerant quantum devices.

2 Method and results

2.1 Overview and input assumptions

To overcome the exponential scaling of the memory cost of predicting generic quantum dynamics by means of NG-RC, we turn to our quantum algorithm for NG-RC. At the heart of NG-RC are the creation of nonlinear feature vectors via tensoring, Eq. 4, and the training to find the optimal (complex-valued) weight matrix W via linear least-square regularization. The predictions of future quantum states are then the results of applying the weight matrix to feature vectors.

Our main quantization techniques are block encoding and the quantum singular value transform (QSVT) (Gilyén et al. 2019; Martyn et al. 2021), which allow us to efficiently implement a series of linear-algebraic operations involving the input and target states that block-encodes W into a unitary gate \(U_W\) as a quantum circuit (middle of Fig. 1). This constitutes the training phase. \(U_W\) is then fed into the circuit of the prediction phase to generate the output state vectors (bottom of Fig. 1).

The input quantum data are fed into the circuits in the form of a coherent superposition of states at different time steps in the dynamics, for example, \(|{o_k}\rangle =(|{s_{k}}\rangle |{0}\rangle +|{s_{k-1}}\rangle |{1}\rangle )/\sqrt{2}\) in the training phase. We formalize this input assumption as having the oracles

$$\begin{aligned} \mathcal {O}_0: |{0}\rangle ^{\otimes d}|{k}\rangle&\mapsto |{s_k}\rangle |{k}\rangle ,\end{aligned}$$
(5a)
$$\begin{aligned} \mathcal {O}_{-1}: |{0}\rangle ^{\otimes d}|{k}\rangle&\mapsto |{s_{k-1}}\rangle |{k}\rangle , \end{aligned}$$
(5b)
$$\begin{aligned} \mathcal {O}_{\tau }: |{0}\rangle ^{\otimes d}|{k}\rangle&\mapsto |{s_{k+\tau }}\rangle |{k}\rangle , \end{aligned}$$
(5c)
$$\begin{aligned} \widetilde{\mathcal {O}}_0: |{0}\rangle ^{\otimes d}|{\tilde{k}}\rangle&\mapsto |{\tilde{s}_{\tilde{k}}}\rangle |{\tilde{k}}\rangle , \end{aligned}$$
(5d)
$$\begin{aligned} \widetilde{\mathcal {O}}_{-1}: |{0}\rangle ^{\otimes d}|{\tilde{k}}\rangle&\mapsto |{\tilde{s}_{\tilde{k}-1}}\rangle |{\tilde{k}}\rangle , \end{aligned}$$
(5e)

where \(|{s_k}\rangle \) and \(|{\tilde{s}_{\tilde{k}}}\rangle \) are d-qubit input data in the training phase and in the prediction phase, respectively, and \(k=0,\ldots ,T-1\) and \(\tilde{k}=0,\ldots ,\widetilde{T}-1\).

To express the total query complexity of the algorithm, we denote the number of calls to the oracle \(\mathcal {O}_i\) (resp. \(\widetilde{\mathcal {O}}_i\)) by \(T_{\mathcal {O}}\) (resp. \(T_{\widetilde{\mathcal {O}}}\)). For convenience, the states created by the oracles are assumed to be error-free, but the total query complexity will take into account errors that occur in the block encoding processes of various matrices.

In practice, the time complexity of the algorithm would depend on the way the oracles are instantiated. For example, we may have an access to the controlled version of a unitary U that generates each data point, a common assumption in the block-encoding literature (Martyn et al. 2021). In this case, even though we could obtain the quantum state in the future \(\tau \) time step by repeatedly applying U, the same state (up to an error) can be generated by a single round of the QNG-RC if it was trained for the time skip of duration \(\tau \).

In general, constructing a controlled-U operation is not possible without further knowledge about U. A weaker but more natural assumption is that we are only handed each data point \(|{s_k}\rangle \), in which case the uniform superposition can still be created by consuming multiple copies of each data point (Kimmel et al. 2017). However, the complexity of the algorithm in Kimmel et al. (2017) depends on the unknown overlap between the data points we wish to superpose.

2.2 Block encoding

We employ the block-encoding technique to construct the non-linear feature matrix and perform the Tikhonov regularization. Specifically, a relevant matrix A is embedded as a submatrix of a unitary gate U such that

$$\begin{aligned} \left( \langle {0}|^{\otimes a}\otimes I^{\otimes s}\right) U \left( |{0}\rangle ^{\otimes a}\otimes I^{\otimes s}\right) = A/\alpha \text { for } U = \left( \begin{array}{ll} A& * \\ *& *\end{array}\right) , \end{aligned}$$
(6)

where \(|{0}\rangle ^{\otimes a}\) is the fiducial state of an a-qubit ancillary system, and \(\alpha \ge \Vert {A}\Vert \) due to the unitary constraint. Definition 1 is a concrete definition of the block-encoding.

Definition 1

Suppose that A is an s-qubit operator, the parameters \(\alpha ,\epsilon \in \mathbb {R}\), and \(a\in \mathbb {N}\). The \((\alpha ,a,\epsilon )\)-block-encoding of a \((s+a)\)-qubit unitary operator U can be defined if

$$\begin{aligned} \Vert {A-\alpha \left( \langle {0}|^{\otimes a}\otimes I\right) U\left( |{0}\rangle ^{\otimes a}\otimes I\right) }\Vert \le \epsilon . \end{aligned}$$

Once the block-encoded matrix is in place, QSVT allows us to construct a degree-q polynomial approximation of essentially any well-behaved function of the singular values of A, using the number of gates U and controlled operations that are efficient in q (Gilyén et al. 2019; Martyn et al. 2021). This approach enables straightforward creation of matrix polynomials and, in particular, creation of the Moore-Penrose pseudo-inverse by inverting the singular values of A. Note that the Tikhonov regularization in the block-encoded form is also presented as a subroutine in Chakraborty et al. (2023). Our approach here, outlined in Lemma S7 in the supplementary material, is slightly distinct from theirs.

2.3 Training phase

According to the NG-RC procedure, by utilizing regularized linear optimization, we obtain an optimal weight matrix through the Tikhonov-regularized pseudo-inverse of the feature matrix X (cf. Eq. 1 with complex elements). Therefore, we begin by efficiently constructing the feature matrix X. Assuming the existence of oracles \(\mathcal {O}_0\) and \(\mathcal {O}_{-1}\), we can generate the linear component of the feature vector, denoted as \(|{o_k}\rangle \), by introducing an additional qubit to entangle with these oracles. This process is depicted by the operator \(U^{\text {lin}}\), shown as a green box in Fig. 1, which maps the state \(|{0}\rangle ^{\otimes d+1}|{k}\rangle \) to \(|{o_k}\rangle |{k}\rangle \). To incorporate the nonlinear component \(|{o_k}\rangle \otimes |{o_k}\rangle \) into the feature vector, we must apply the operator \(U^{\text {lin}}\) twice due to the constraints imposed by the no-cloning theorem. The resulting feature vector \(|{x_k}\rangle \) Eq. 4 is represented by the quantum circuit \(U^f\), which maps the state \(|{0}\rangle ^{\otimes 2d+3}|{k}\rangle \) to \(|{x_k}\rangle |{k}\rangle \), as also illustrated in Fig. 1 (middle, orange box). Note that, with this encoding method, the dimension of \(|{x_k}\rangle \) is \(8D^2\), wherein \(4D^2+2D\) dimensions constitute the feature vector \(\varvec{x}_k\) analogous to the case of classical NG-RC Eq. 4, and the remaining dimensions are for zero padding. We can coherently construct the feature matrix \(X=\left( |{x_0}\rangle ,\cdots ,|{x_{T-1}}\rangle \right) \in \mathbb {C}^{8D^2\times T}\) with the block-encoding technology using quantum gates

$$\begin{aligned} \left( I^{\otimes \max (0,2d+3-t)}\otimes H^{\otimes t}\otimes I^{\otimes \max (2d+3,t)}\right) \cdot \text {SWAP}\cdot U^f, \end{aligned}$$

which is the \((\sqrt{T},\max (2d+3,t),0)\)-block-encoding of the feature matrix X, illustrated in Fig. 1 (middle). Note that the dimension of X here is \(\max (T,8D^2)\times \max (T,8D^2)\), which includes another zero padding to convert the submatrix X in the block-encoding unitary into a square submatrix (Gilyén et al. 2019).

Due to the requirement for matrix multiplication involving Y and the pseudo-inverse of X in Eq. 1, the size of submatrix Y in the block-encoding unitary must align with the dimension of the submatrix X. Utilizing block encoding techniques, the quantum circuit

$$\begin{aligned} \left( I^{\otimes \max (0,2d+3-t)}\otimes H^{\otimes t}\otimes I^{\otimes \max (2d+3,t)}\right) \cdot \textsc {SWAP}\cdot \mathcal {O}_\tau , \end{aligned}$$

along with pre-amplification proposed in Chakraborty et al. (2019), can effectively create the \((\sqrt{2}\Vert Y\Vert ,\max (2d+3,t)+1,\delta _Y)\)-block encoding of Y, where \(\delta _Y\) is the error of the pre-amplification.

By virtue of Lemma S7 in the supplementary material, the optimal regularized weight matrix Eq. 1 in the block-encoded form can be constructed with the resources stated below in Theorem 2.

Theorem 2

Let \(\delta _W\in (0,1]\). Suppose that we have the \((\sqrt{T},\max (2d+3,t),0)\)-block-encoding of the feature matrix X and the \((\sqrt{2}\Vert Y\Vert ,\max (2d+3,t)+1,\delta _Y)\)-block encoding of training target matrix Y. Let \(\kappa _X\) be the condition number of X, and

$$\begin{aligned} \kappa = \kappa _X\sqrt{\frac{\Vert {X}\Vert ^2+\lambda }{\Vert {X}\Vert ^2+\lambda \kappa _X^2}}. \end{aligned}$$
(7)

where \(\lambda \ge 0\) is the regularization parameter (Chakraborty et al. 2023). Also, define

$$\begin{aligned} w=2\max (2d+3,t)+2 \end{aligned}$$
(8)

to be the number of ancilla qubits used in the block encoding of W. In a condition such that \(\delta _Y\le \delta _W/(4\kappa )\), then we can construct the unitary operator \(U_W\) as a \(\left( \frac{2\sqrt{2}\kappa \Vert Y\Vert }{\Vert {X}\Vert +\sqrt{\lambda }},w,\delta _W\right) \)-block-encoding of the weight matrix W in

$$\begin{aligned} T_W=O\left( \left( \frac{\kappa }{\Vert {X}\Vert +\sqrt{\lambda }}+\frac{1}{\Vert {Y}\Vert }\right) \log \left( \frac{\kappa \Vert {Y}\Vert }{\delta _W}\right) T_\mathcal {O}\sqrt{T}\right) \end{aligned}$$
(9)

queries, where \(T_\mathcal {O}\) is the number of calls to the oracles.

In the case \(t>2d+3\), \(U_W\) would be a \(\Big (\frac{2\sqrt{2}\kappa \Vert Y\Vert }{\Vert {X}\Vert +\sqrt{\lambda }},2t+2,\) \(\delta _W\Big )\)-block encoding of the weight matrix W. The query complexity Eq. 9 can be expressed more directly in terms of T and D using inequality of operator norm. Since each column in both X and Y is unit vector, implying that we have an inequality \(\frac{1}{\sqrt{T}}\le \Vert {X}\Vert \le D\sqrt{T}\) and \(\frac{1}{\sqrt{T}}\le \Vert {Y}\Vert \le \sqrt{DT}\), respectively. Matrix norm inequalities are detailed in Section S2 B of the supplementary material. By assuming \(\lambda =O(1)\), the time complexity is expressed to

$$\begin{aligned} T_W = O\left( \kappa T\log \left( \frac{\kappa DT}{\delta _W}\right) T_\mathcal {O}\right) , \end{aligned}$$
(10)

Note that the time complexity here is difficult to compare to that of the classical NG-RC because of the unknown gate complexity of implementing \(\mathcal {O}_i\). We can only say that the number of calls to \(\mathcal {O}_i\) is of the order of \(O\left( \kappa T\log \left( \frac{\kappa DT}{\delta _W}\right) \right) \).

Fig. 2
figure 2

The performances of the NG-RC in predicting unseen future states of a four-qubit transverse-field Ising model in the disordered phase with \(J=0.5\) and \(h=5\) across \(\widetilde{T} = 4 \times 10^4\) future time steps employing two different approaches: the skipping-ahead method with a time skip of \(\tau =10^6\) steps (the map \(f_\tau \) represents the trained NG-RC for predicting the next \(\tau \) step from the current step), and the conventional approach of iteratively predicting each successive time step (\(g^{(\tau )}\) is the composition of \(\tau \) successive g obtained from the trained NG-RC \(f_{\tau =1}\). \(T=2\times 10^4\) steps of the time evolution with the same step size of \(\Delta t=1/(200E_{\text {max}})\), where \(E_{\max }\) is the highest eigen-energy of the system, are used to train the NG-RC in both cases. The training is optimized with the optimal regularization parameter \(\lambda = 1.0\times 10^{-3}\) in the iterative prediction and with \(\lambda = 0\) in the skipping-ahead prediction. Our benchmark targets are the real and the imaginary parts of all amplitudes in the computational basis ((a) and (b)). The comparisons of the expectation values \(\langle {X_0}\rangle \), \(\langle {X_0X_1}\rangle \), and the fidelity \(F(\tilde{s}_t,s_t):=|\langle {\tilde{s}_t}|{s_t}\rangle |\) between the target and the predicted states are shown in (c) and (d)

Full size image

2.4 Prediction phase

Now that we have the trained weight matrix in the block-encoded form, \(U_W\), we can build a circuit to predict future quantum states by the skipping-ahead method. As shown in Fig. 1, the quantum circuit for the prediction phase of the QNG-RC has the same structure as that for the training phase, but with input past quantum states incorporated via the oracles \(\widetilde{\mathcal {O}}_0\) and \(\widetilde{\mathcal {O}}_{-1}\) defined in Eqs. 5d and 5e, respectively. The output of the circuit is the desired prediction vector \(|{\tilde{s}_{\tilde{k}+\tau }}\rangle =W|{\tilde{x}_{\tilde{k}}}\rangle /\Vert {W|{\tilde{x}_{\tilde{k}}}\rangle \Vert }_2\) up to an error \(\delta \in (0,1)\) due to an inaccuracy in the application of the weight matrix W. More specifically, we have an inequality

$$\begin{aligned} \delta _W\le \frac{\delta \Vert {W}\Vert }{4\kappa _W}, \end{aligned}$$
(11)

where \(\delta _W\) is the error in Theorem 2, and \(\kappa _W\) is the condition number of W. Employing the pre-amplification technique for block encoding applied on a quantum state in Ref. Chakraborty et al. (2023), the number of queries to produce such a \(\delta \)-close prediction vector is

$$\begin{aligned} O\left( \frac{\kappa \kappa _W}{\Vert {X}\Vert +\sqrt{\lambda }}\frac{\Vert {Y}\Vert }{\Vert {W}\Vert }\log \left( \frac{\kappa _W}{\delta }\right) T_W+\kappa _WT_{\widetilde{\mathcal {O}}}\right) , \end{aligned}$$

where \(\kappa \) is defined in Eq. 7. What is more, the operator norm inequality allows us to rewrite the query complexity as \(O\left( {\kappa \kappa _W}\log \left( \kappa _W/\delta \right) T_W+\kappa _WT_{\widetilde{\mathcal {O}}}\right) \). The technical details regarding the pre-amplification and the operator norm inequality are given in Lemma S4 and Section S2 B in the supplementary material.

Number of ancilla qubits.–The overall circuit for the prediction represents the unitary operator

$$\begin{aligned} \tilde{U}: |{0}\rangle ^{\otimes w'}|{0}\rangle ^{\otimes 2d+3}|{\tilde{k}}\rangle \mapsto |{0}\rangle ^{\otimes (w'+d+3)}|{\tilde{s}_{\tilde{k}+\tau }}\rangle |{\tilde{k}}\rangle , \end{aligned}$$
(12)

where the register \(|{\tilde{k}}\rangle \) contains \(\max \left( \tilde{t},2d+3\right) \) qubits. The block encoding requires in total \(w'=w+\max (0,t-2d-3)\) qubits. The w given in Theorem 2 is the number of index qubits in the block-encoding (BE) register, and \(\max (0,t-2d-3)\) is the number of qubits used in the zero padding in the register RD, see Fig. 1 and the corresponding caption. In the case \(t>2d+3\), the number of ancilla qubits is \(w'=3t-2d-1\) and \(w'=4d+8\) for \(t<2d+3\).

3 Conclusion and outlooks

Drawing inspiration from the paradigm of NG-RC, we develop a novel, end-to-end quantum algorithm for predicting skipped-ahead many-body quantum dynamics. The algorithm requires a time-series of quantum data with no assumption on the unitary generating the dynamics, i.e., the forecasting is model-free. The algorithm employs block encoding to efficiently handle matrix algebra subroutines, including matrix multiplication, inversion, and regularized linear optimization, which are manipulated by the QSVT. In addition to providing a quantum computational speedup and avoiding exponential resource consumption in storing many-body quantum states classically, our algorithm coherently processes and generates quantum data, thereby circumventing classical-quantum data conversion problems during encoding and readout procedures. Since the main advantage of the QNG-RC is the coherent output of far-future quantum states, it is worth investigating how the QNG-RC can be interfaced with or used alongside other algorithms. One potential application could be using the QNG-RC to output future states within a thermalized quantum ensemble, thereby functioning as a sampling algorithm. However, the specific dynamics and initial states for which the QNG-RC is most effective in such scenarios remain an open question and warrant further studies.

Both the NG-RC and the QNG-RC algorithms are heuristic in the sense that the convergence of the predictions to the true state is not guaranteed. Therefore, determining the effectiveness of this method for larger quantum systems and developing benchmarks when classical simulation is no longer feasible remains an outstanding challenge. In contrast, it is worth noting that our algorithm is not heuristic during the training stage. The “trained” optimal weight matrix is derived as a closed-form solution of the regularized least-squares optimization problem, rather than through gradient descent. As a result, the training is not susceptible to issues such as the barren plateau phenomenon that often affects the training of variational quantum circuits (McClean et al. 2018).

In the original NG-RC framework, a feature vector \(\varvec{x}_i\) can be constructed with a degree-p monomial of an m-time-delay data, see Eqs. 2–3. In our quantum circuit design, though we assumed \(m=2\) for simplicity, the quantum operator \(U^\text {lin}\) can be slightly modified to incorporate the oracles \(\{\mathcal {O}_{-j\Delta }\}_{j=0}^{m-1}\) for an m-time-delayed \(|{o_k}\rangle \). A quantum circuit that encodes the nonlinear map of Fig. 1 would then involve the addition of multiple \(U^\text {lin}\) and more qubits to accommodate the tensor product \((\varvec{o}_i)^{\otimes p}\). The construction of the corresponding quantum circuits of this more general form of non-linearity is elaborated in Section S3 of the supplementary material.

Lastly, the conventional application of the NG-RC to make iterative predictions, whose results are shown in Fig. 2b, can also be implemented on a quantum circuit. However, the number of gates and the error scales exponentially with the number of iterations. Nevertheless, we discuss a quantum circuit for this approach in Section S4 of the supplementary material, along with the analysis of the error propagation.