Estimating many properties of a quantum state via quantum reservoir processing

Estimating properties of a quantum state is an indispensable task in various applications of quantum information processing. To predict properties in the post-processing stage, it is inherent to first perceive the quantum state with a measurement protocol and store the information acquired. In this work, we propose a general framework for constructing classical approximations of arbitrary quantum states with quantum reservoirs. A key advantage of our method is that only a single local measurement setting is required for estimating arbitrary properties, while most of the previous methods need exponentially increasing number of measurement settings. To estimate $M$ properties simultaneously, the size of the classical approximation scales as $\ln M$ . Moreover, this estimation scheme is extendable to higher-dimensional systems and hybrid systems with non-identical local dimensions, which makes it exceptionally generic. We support our theoretical findings with extensive numerical simulations.


I. INTRODUCTION
Estimating properties of a quantum state plays a central role in the implementation of various quantum technologies, such as quantum computing, quantum communication and quantum sensing.This highlights that extracting information from a quantum system to a classical machine lies at the heart of quantum physics [1].The prominent technique for this task, quantum tomography, studies the reconstruction methods of density matrix, which captures all the information of a quantum system.However, the curse of dimensionality has emerged with the advent of the noisy intermediate-scale quantum (NISQ) era [2], which renders it infeasible to obtain a complete description of quantum systems with a large number of constituents.Moreover, a full description is often superfluous in tasks where only key properties are relevant.As a consequence, the concept of shadow tomography is proposed to focus on predicting certain properties of a quantum system [3].
A particularly important progress in the study of shadow tomography is the advancement of randomized measurements [4][5][6], the virtue of which is highlighted as "Measure first, ask questions later" [7].The randomized measurement protocols proposed by Huang, Kueng and Preskill construct approximate representations of the quantum system, namely classical shadows, via Pauli group and Clifford group measurements [5].The single-snapshot variance upper bound of classical shadows is determined by the socalled shadow norm, which is optimal in the worstcase scenario.Nevertheless, variants of randomized measurements can achieve a better performance in specific cases [8][9][10].For instance, the Hamiltoniandriven shadow tomography achieves a higher efficiency in predicting diagonal observables [10].In addition, the readout noise present in the measurement protocol can be further suppressed by constructing classical shadows with positive operator-valued measures (POVMs) [11][12][13].
The classical shadows are highly efficient in the estimation of various properties in the post-processing phase, the benefits of which extend to entanglement detection [14], characterization of topological order [15], machine learning for many-body problems [16], etc.However, the randomized measurements protocols pose a challenge in experiments due to the need for exponentially increasing measurement settings to achieve an arbitrary accuracy.Hence, various techniques are introduced to tackle this problem [17][18][19][20].Moreover, the theoretical results are based on the fact that multi-qubit Clifford groups are unitary 3-designs [21], which is not the case for arbitrary qudit systems.The generalization of these results to higher-dimensional systems typically require complex unitary ensembles that are relatively hard to implement [22,23].Therefore, a general method for direct estimation with a single measurement setting is highly desirable.
Recently, quantum neural networks [24][25][26] are widely studied as promising artificial neural networks due to their enhanced information feature space supported by the exponentially large Hilbert space [27,28].Unlike traditional computing frameworks, neural networks learn to perform complex tasks based on training rather than predefined algorithms or strategies [29].With the capacity to produce data that displays atypical statistical patterns, quantum neural networks have the potential to outperform their classical counterparts [25].However, training a quantum neural network can be equally hard [30].Indeed, it has been shown that training of quantum neural networks could be exceptionally difficult owing to the barren plateaus or far local minima in the training landscapes [31][32][33][34].This is the reason that quantum neural networks are often limited to shallow circuit depths or small number of qubits.A trending line of research that circumvents this issue is quantum reservoir processing (QRP) or reservoir computing [35][36][37][38][39], which studies the quantum analogy of recurrent networks.Note that in this context, 'reservoir' refers to a type of neural network.In QRP, training is completely moved out of the main network to a single output layer, such that the training becomes a linear regression eliminating the possibility of producing barren plateaus or local minima [35].Such a quantum neural network retains its quantum enhanced feature space while being trainable via a fast and easy mechanism.
In this work we present a shadow estimation scheme on the QRP platform, which overcomes the obstacles faced by randomized measurement protocols by harnessing the richness of QRP.While QRP serves as a quantum analog of reservoir computing, our focus in this study is limited to its role as a quantum state processing platform that extracts properties of the input quantum system.Therefore, the training at the output layer maintains as a straightforward approach by utilizing linear inversion.A scheme of minimal quantum hardware comprising pair-wise connected quantum nodes is developed to estimate many properties of a quantum state.We highlight that only two-node training data is required, which captures the pair-wise interacting reservoir dynamics.As major advantages, our scheme requires single-node measurements, only in a single setting, and a linear reservoir size 2n with respect to the number of constituents n of the input state.All of these are particularly favorable for actual physical implementations.Furthermore, we establish rigorous performance guarantee by adopting the mindset of shadow estimation.According to Born's rule, one measurement of a quantum state is analogous to sampling a probability distribution once.Thus, learning properties of a quantum state involves measuring identical and independently distributed (i.i.d.) samples of the quantum state a certain number of times.To estimate M observables of the state within an additive error ϵ and with constant confidence, the number of i.i.d.input samples consumed scales as O F res ln M/ϵ 2 .The factor F res represents the variance upper bound of the single sample estimator, which depends solely on the observables and the reservoir dynamics, and its magnitude is comparable to that of the shadow norm.As a direct consequence of the pair-wise reservoir dynamics, F res for a k-local tensor product observable is the product of that for each single-qubit observable.We support the theoretical results with extensive numerical simulations.

II. QUANTUM RESERVOIR PROPERTY ESTIMATION
In this section we introduce the scheme of quantum reservoir property estimation (QRPE), which evaluates physical properties of input quantum states based on a quantum reservoir processing device.Our goal parallels that of shadow estimation: to devise a resource-efficient classical representation of complex quantum states that permits access to their properties through subsequent classical processing.

A. Physical setting
As a proof-of-principle, we first consider a dynamical estimation device based on pair-wise connected interacting qubits, as shown in Fig. 1, which is also known as quantum registers [40,41].Later in Sec.IV we will demonstrate that our scheme works also for continuous variable reservoir systems evolving under open quantum dynamics.The connections between quantum nodes are obtained with transverse exchange interactions [42] and each qubit is excited with an onsite driving field.For n-qubit input states, there are n pairs of reservoir nodes.The corresponding Hamiltonian of the device is given by The operators ς x,y,z i represent the Pauli operators on the i-th quantum node, which has a compatible di-mension with the context.The parameter J represents the strength of the pair-wise transverse exchange interaction between the reservoir nodes.The parameters P 1,2 and E 1,2 represent driving field strength and onsite energy respectively.Such a device can be readily realized with superconducting qubits, where the exchange interaction can be realized via a cavity quantum bus [42].Moreover, the form of Ĥ is of a quantum spin Hamiltonian which can be realized in a variety of platforms such as NMR [43,44], quantum dots [45], and trapped ions [46].

B. Measurement protocol
The quantum state of interest, which we denote by σ, is injected into the quantum reservoir via an invertible map.For simplicity, we consider the swap operations.For a pair of interacting nodes ⟨2i − 1, 2i⟩, the (2i−1)th reservoir node is connected to the input system.Thus, the initial state of the network at time t = 0 is given by the density matrix: where the suffix 'rest' indicates the network nodes other than the ones connected to the input qubits via the swap gates.The initial reservoir state evolves in time as, where ρ(t) is the density operator at time t and Û (t) = exp(−it Ĥ/ℏ) is the evolution operator.After a sufficient time evolution, we perform local Pauli-Z measurements on the reservoir nodes (qubits).For each node there are two readouts, +1 and −1, represented by the projectors onto the positive and negative eigen-subspace of Pauli-Z operator {ς z i } respectively.The final readouts are provided by a set of commuting readout operators where N node is the number of reservoir nodes.ĈS is an element of the set {ô i }, which is related to the Pauli-Z operators as so it represents the configuration of measurement outcomes that only nodes in the set S result in −1 for local Pauli-Z measurements.Hence, the readout for the i-th input is recorded by a bit string s i ∈ {0, 1} n .The total number of readout operators in {ô i } is given by where N node k is the combinatorial number.
Considering that {ô i } forms an orthogonal measurement basis with a total number of 2 N node elements, a reservoir with a minimal of N node = 2 log 2 d nodes is required to estimate arbitrary quantum properties for input states supported on a d-dimensional Hilbert space H d .This observation agrees with our proposal of pair-wise connected reservoir networks.Moreover, the commutativity of the chosen readout operators leads to quantum resource effectiveness.This is in sharp contrast to the traditional quantum reservoir computing schemes where either the required size of the quantum reservoir or the temporal resolution in the measurement tend to be exponentially large (∼ d 2 ).In either situations, these traditional schemes are exponentially quantum resource consuming.

C. Training
Here we introduce a training process that captures the internal maps of QRP.In the context of QRP, the essential condition is that we have reproducible reservoir dynamics.However, the reservoir could be largely a black box, especially when we consider open quantum dynamics.Thus, a training process at a single output layer is required.The basic property of a quantum state σ is the linear function in the form of Tr(Oσ), where O is a linear operator that is compatible with σ.
For instance, the expectation value of an arbitrary observable O takes on this form.To estimate the linear function Tr(Oσ), the QRPE scheme essentially maps the input state σ to a vector of probabilities for observing each readout operator and the target observable O to a vector of weights satisfying We note that similar maps have also been studied in the context of analog quantum simulation [19,20].
Here each readout in experiments requires only linear classical storage with respect to the system size, owing to the tensor product structure in Eq. ( 3).The relation given by Eq. ( 6) is achieved by sampling from i.i.d.copies of the n-qubit state σ and processing the reservoir readouts with statistical methods, as addressed in Sec.II D, while Eq. ( 7) is achieved by a training process described below.
For training, we require a one-time estimation of a known set of training states {|φ k ⟩}.For the simplest qubit-node reservoir system as described in this section, there is no setback in having a complete characterization over the reservoir dynamics.However, when the reservoir dynamics is much more complex as discussed in Sec.IV, the training process described here can still capture the reservoir maps with a resource cost depending only on the dimension of input quantum system.We note that the training process is similar to a quantum process tomography over a small quantum system.
Until this step, the whole procedure is completely independent of the property to be estimated.The knowledge of O is only required at the post processing level, where we set the target output Let where Y out is a row vector with elements Y out k , k = 1, 2, . . ., N train , and Y tar is defined similarly.For typical quantum reservoirs, a total of d 2 training states are needed to estimate arbitrary properties for an input state supported on H d .The set of training states {|φ k ⟩⟨φ k |} is a set of d 2 vectors which spans a d 2 -dimensional Hilbert space, i.e., forms an informationally complete POVM.We choose the training states as where where Then the expectation of reservoir readout for an input state σ is as is explained in Appendix.A. If X t is full rank, there exists a weight vector W that minimizes E train , i.e., which can be written as W = ⟨⟨O|T −1 .
Leveraging pair-wise reservoir dynamics in Eq. ( 1) lifts the burden on training, given that the training data inherits a tensor product structure.Once the training data of a single pair of interacting reservoir nodes is collected as X p , then we have where n is the number of node pairs.Also, where T p = X p M −1 p .Thus, the training task is effectively reduced to that of a two-node reservoir, and the full-rank requirement of X t is correspondingly reduced to that of X p .Moreover, the vector of weights only necessitates polynomial storage if the property O can be decomposed into a finite sum of tensor products.See Appendix.A for more details.

D. Reservoir estimator
In this section we introduce the reservoir estimators for linear functions.With the training data at our disposal, we could analyze the sample efficiency of the reservoir estimators.
Suppose the observed readout operator for the i-th input copy is ôj , then the so-called single snapshot X i is a vector where the j-th element is 1 and the other elements are 0.For a total of N sample input copies, one obtains a set of snapshots where X is a random variable that conforms to the probability distribution behind X i .Eq. ( 8) indicates that Ω is an unbiased estimator for Tr(Oσ).Hence in data processing, we could apply the median of means (MoM) method to neutralize the effect of outliers [5,8].After processing N sample = KN input copies, we divide the snapshots into K equally sized subsets {X v i | i = 1, 2, . . ., N }, and compute the mean value of the single-snapshot estimators for each subset.The corresponding estimators are Then, the MoM estimator is given by With this, we have Theorem 1.The number of quantum state inputs needed for estimating a set of M properties {O i | i = 1, 2, . . ., M } to precision ϵ and confidence level 1 − δ scales in The factor F i res = ||B i || ∞ is the variance upper bound of the single-snapshot estimator of O i maximized over the possible quantum state inputs, where || • || ∞ represents the spectral norm, and B i is defined by Proof.This efficiency scaling results from the median of means method, and we consider the worst-case scenario by maximizing F i res over the set of observables.A more detailed proof is included in Appendix B.
The variance of a single-snapshot estimator for O is invariant for properties Thus, one could use only the traceless part of O to compute the worst-case variance upper bound.It is interesting to note that the MoM estimator won't have a visible significance in some tested cases [47,48], where it could be replaced with the sample mean estimator.The performance factor F res allows for an investigation of the reservoir parameters in Eq. (1).To achieve tomographic completeness, it is essential to have sufficiently large evolution time t and hopping strength J, ensuring effective information scrambling.We also find that when both E 1 and E 2 equal 0 or both P 1 and P 2 equal 0, the reservoir map will be tomographically incomplete, which is closely related to the measurement setting.If we replace ς z with ς x or ς y in Eq. ( 3), then the reservoir map can be tomographically complete when both P 1 and P 2 equal 0. However, the reservoir settings with a good performance seem to happen without a pattern in numerical experiments.The performance of one reservoir setting at different time and random reservoir settings is illustrated in Fig. 8 and Fig. 9 respectively in Appendix.B. We also find that increasing more nodes or non-local interactions is unlikely to improve the performance.On the contrary, in numerical experiments of reservoirs with random nearest-neighbor interaction or fully-connected interaction, the average performance is worse than that of the pair-wise interacting reservoir by several orders of magnitude.
To further analyze the sample efficiency scaling of the QRPE scheme, we have: Theorem 2. For a k-local tensor product observable, e.g. the worst-case variance upper bound ||B|| ∞ is the product of that of each local observables O i , i.e., where B i satisfies Proof.This theorem results from the pair-wise reservoir dynamics.See Appendix B for the details.
Theorem 2 indicates that for k-local tensor product property estimation, if a reservoir setting works well in the single-qubit case, then it also works well in the multi-qubit case.Thus, our approach for evaluating the reservoir parameters J, P 1,2 , and E 1,2 is based on the performance in the single-qubit state overlap estimation task.To see why overlap estimation reflects the overall performance of observable estimation, note that the single-snapshot estimator's variance of an arbitrary property O ′ satisfying O ′ = c 1 σ + c 2 1 1 equals c 2 1 times that of σ, where c 1,2 are real numbers and σ is a density matrix.We find several settings that lead to small average variance upper bound of random target states, and choose one of them as the reservoir setting used in this manuscript.
In the task of pure state fidelity estimation, the efficiency of the current reservoir estimator outperforms the random Pauli measurements only in a fraction of target pure states.We note that the optimal sample efficiency of shadow estimation with local measurements is achieved with random Pauli measurements [5,11].Also, the ratio of the average variance upper bound of shadow estimation and that of the QRPE scheme is slightly larger than unity and almost invariant with the system size.While the average efficiency of the current setting only marginally differ from that of the shadow estimation, the number of settings for the standard shadow estimation is exponentially large compared to that of the present scheme.These results are shown in Fig. 2. While one could obtain a scaling independent of the system size with global Clifford group measurements, it would require to implement complex Clifford compiling [50,51].The reservoir estimators for nonlinear functions are based on U-statistics [5,52], which is a generalization of the sample mean estimator.It is worth noting that the reservoir snapshots are ready to be used to estimate future properties of the input state σ.

III. APPLICATIONS
Here we provide a wide range of applications of the QRPE scheme.The reservoir initialization and evolution time are fixed for all applications.See Fig. 2 for the details of reservoir settings.

A. Fidelity estimation
Quantum fidelity is a widely used distance measure for quantum states, which is a linear function of the given state when the target state is a pure state [53,54].In Fig. 2, we present an illustrative example of how the reservoir estimation scheme works in fidelity estimation.Numerical results indicate that for the average worst-case variance upper bound for pure target states randomly generated from Haar measure, the ratio between QRPE and the random Pauli measurement is close to 1, while there exists an exponential improvement in the number of measurement settings.

B. Entanglement detection
Entanglement is an indispensable resource in tasks ranging from quantum computation to quantum communication [55].Whether one can claim its existence for a given system has aroused both theoretical and experimental interests [56][57][58][59][60].However, the difficulty in describing the convex space of separable states poses a trade-off relation between the detection ability and effectiveness of entanglement criteria [61].With its capability of estimating multiple observables simultaneously, the QRPE scheme is a natural fit for the task of entanglement detection with linear entanglement criteria.
We illustrate the QRPE scheme on the detection of an intriguing and more targeted phenomenon, the entanglement sudden death (ESD) [62,63].Consider a three-qubit GHZ-type state [64] where ρ G = (|000⟩ + |111⟩)(⟨000| + ⟨111|)/2.The dephasing channel is defined by the Kraus operators Set p = 1 − exp (−κt), then dephasing of the initial state ρ reflects on a factor exp (−κt) that times the off-diagonal elements.In Fig. 3, we use the following optimal linear entanglement witnesses for detecting genuine multipartite entanglement (GME) and any multipartite entanglement (ME) respectively [65]: where ρ G− = (|000⟩ − |111⟩)(⟨000| − ⟨111|)/2.We normalize the spectral norm of the two witnesses, and estimate both of them simultaneously.With a total of 6000 snapshots for each {q, κt}, the maximal difference from the true value is 0.13.We see that under a small dephasing noise, both the GME and ME witnesses vanish at a finite time for some initially entangled states, indicating the emergence of ESD.

C. Estimating expectation values of local and global observables
Estimating the values of observables is a fundamental task in various quantum information processing tasks.For local observables, we present the numerical results of estimating the local Pauli observables of randomly chosen four-qubit states.The numerical result is shown in Fig. 4. We observe that the sample complexity required for achieving a given error rate of estimating Pauli observables increases exponentially with the locality, which agrees with Theorem.2.
For global observables, we consider the task of GHZ state fidelity estimation.For an input GHZ state with three, six, nine and twelve qubits, we present the results of numerical experiment as the number of input copies versus the average error in Fig. 5.It can be seen that the sample complexity of identifying GHZ states with a given error rate increases exponentially with the number of qubits, which is similar to the case of local tensor product observables.

D. Estimating nonlinear functions
The reservoir measurements only describe linear functions in Eq. ( 8).Nevertheless, experimentally accessible nonlinear functions are typically measured by es-timating linear functions that can be translated into the QRPE scheme.For instance, the swap trick is widely used in purity estimation.For two copies of the input state σ, the swap operator S obeys Tr(Sσ ⊗ σ) = Tr(σ 2 ).Numerical simulation of the purity estimation of 10000 random one-qubit input states shows that the average variance upper bound is around 2.9.
The Rényi entropy is another important nonlinear factor for characterizing entanglement, which is the logarithm of subsystem purity.The estimator of second Rényi entropy is where S A is the swap operator acting on subsystem A of two copies of ρ.The second Rényi entropy of small subsystems is useful in avoiding weak barren plateaus (WBP) [66].Here we perform WBP diagnosis with the reservoir estimation scheme.For an initial state |0⟩ ⊗14 , each gate sequence of the variational quantum eigensolver (VQE) circuit is composed of random local rotations exp(− i 2 θς), where θ ∈ [π/20, π/20] and ς ∈ {ς x , ς y , ς z }, and nearest-neighbor controlled-Z gates with periodic condition.We detect the emergence of WBP by estimating the second Rényi entropy of a small region consisting of N A qubits via 2000 measurements at each circuit depth.The numerical experiment is shown in Fig. 6.We observe that the shaded region created by 10 independent estimations is centered by the real value of the second Rényi entropy with a small fluctuation, which grows with the size of the subsystem.Within a circuit depth of 100, the QRPE scheme efficiently detects the emergence of WBP by the clear phenomenon that the estimated value of second Rényi entropy approaches the page entropy.

E. Higher-Dimensional and Hybrid Systems
Higher-dimensional systems and hybrid systems with non-identical local dimensions are of fundamental importance as a playground to reveal interesting quantum phenomena, such as quantum steering [67,68] and the demonstration of contextuality and nonlocality trade-off [69,70].Here we demonstrate the flexibility of our scheme beyond the qubit systems.
An intuitive generalization from the qubit case when the input state exists in a higher-dimensional or hybrid system is to design a quantum reservoir that has a compatible dimension.For example, if the input state consists of a qubit and a qutrit, then a quantum reservoir with two interacting pairs of qubits and qutrits can readily estimate properties of it.The bosonic Hamiltonian for a pair of interacting reservoir nodes where the operators â represent lowering operators of the quantum nodes (qudits).The parameters J, P 1,2 , E 1,2 and α 1,2 represent hopping, onsite driving field, onsite energy and nonlinear strength respectively.Similar to the qubit case, the reservoir map will be tomographically incomplete when P 1 and P 2 both equal 0 or E 1 and E 2 both equal 0. For each pair of reservoir nodes ⟨i, j⟩, the i-th node is connected to the input system, and the number of pairs equals the number of constituents of the input state.We consider the unitary time evolution governed by the quantum Liouville equation.
For qudit systems with local dimension d, a superoperator basis consists of the generalized Gell-Mann matrices and the normalized identity [71].The readout operators are constructed in the same way with that of the qubit system, the only difference is that there are d projection operators instead of two, corresponding to the d possible population numbers.Due to the tensor product structure of the measurements and training states, the reservoir estimation scheme can be naturally extended to hybrid systems with non-identical local dimensions, such as qubit-qutrit systems.The pair-wise reservoir setting for a pair of qudit nodes is the same with that of qudit reservoirs.
For application we apply virtual distillation (VD) [72][73][74] to estimate the fidelity of a noisy state ρ ϵ , where The estimator for ρ m ϵ is chosen as where ρsi is the single-snapshot estimator of the noisy state ρ ϵ , P (N, m) is the m-permutations of N , * denotes the summation over all distinct subscripts, i.e., s 1 , s 2 , . . ., s m is a m-tuple of indices from the set {1, . . ., N } with distinct entries.The denominator is estimated similarly.The estimation results are illustrated in Fig. 7 for the following maximally entangled qubit-qutrit and two-qutrit pairs, We observe that the virtually distilled states exhibit a significantly improved fidelity as compared to the nondistilled states.Moreover, the statistical fluctuation for the 2 × 3 system is smaller than that of the 3 × 3 system, which is consistent with the performance in single-qudit pure state fidelity estimation.In the given reservoir setting for qutrit systems, the average F res for estimating fidelity of randomly chosen singlequtrit pure state is around 2.76, which is worse than the qubit case.

IV. CONTINUOUS VARIABLE QUANTUM RESERVOIRS
In this section we study a bosonic quantum reservoir that exists in continuous variable (CV) systems, while the input state emitted by a source still lies in discrete variable (DV) systems.The reservoir consists of identical pairs of interacting nodes, as shown in Fig. 1.A single pair of interacting reservoir nodes ⟨k, l⟩ is represented by the Hamiltonian where âi is the field operator of the i-th reservoir node, P represents the strength of the coherent driving field, E i is the onsite energy and α is the onsite interaction strength (Kerr type nonlinearity).We use essentially the same readout and training process as for the DV quantum reservoirs.Here the readout operators are projectors of Fock basis states.For instance, if the input state exists in a qubit system, then the measurements at each reservoir node include the first Fock state |0⟩⟨0| and its complementary operator 1 1−|0⟩⟨0|.
For a single pair of CV reservoir nodes, the input single-qudit system exists in a bosonic mode b1 .The density matrix ρ of the entire system is in both of the quantum reservoir and incident modes, which is governed by the quantum master equation [35,75,76] where is the Lindblad operator acting on the field operator x; γ represents the decay rate; W in represents the input weight and η = (W in ) 2 is set to remove the source photons that have excited the reservoir.On the right side of Eq. ( 35), the first term is the coherent Hamiltonian evolution of the reservoir; the second term represents the decay of reservoir modes; the second line represents the cascade coupling between the input mode b1 and the reservoir mode â1 ; and the last line represents the decay of the input mode.At time t ∈ [t 1 , t 1 + τ ], f (t) = 1 and the input mode b1 is connected to the reservoir node.At other time points f (t) equals zero.We measure the reservoir nodes at time t = t1 + τ + t2.The quantum master equation can be solved with the numerical method provided in Ref. [36].Additionally, Eq. ( 35) can be naturally generalized to the case where the input system has multiple constituents, each of which undergoes the same dynamical process.
In numerical simulation, the truncated Fock space we use has a maximal photon number of 10.Here we set t 1 = 0.3, τ = 1.4,t 2 = 0, γ 1 = 0.53, γ 2 = 0.87, P = 0.95, W in = 0.56, E 1 = 0.94, E 2 = 0.57, J = −0.22 and ω = 1.We choose α i /γ i = 0.15, a ratio consistent with the experiment results [77].The spectral radius [36], which is the largest eigenvalue of the hopping part of the reservoir Hamiltonian, is set to 1.5.With this reservoir setting, the average F res for estimating fidelity of randomly chosen single-qubit pure state is around 9.16.With numerical results, we find that the reservoir map becomes tomographically incomplete when the onsite driving field of the input system is missing.On the contrary, tomographic completeness can still be achieved if we set E i , ω and α equal 0 or neglect the onsite driving field of the reservoir system.Also, we observe that the impact of Kerr nonlinearity on performance in numerical tests is not substantial.Given the infinite dimensionality of the CV reservoir system, it has the capability to process any higher-dimensional or hybrid input states.However, with the same reservoir setting as described above, the average F res for estimating fidelity of randomly chosen single-qutrit pure states is around 490.Thus, similar to the DV quantum reservoir, the performance of CV quantum reservoir settings largely depends on the dimension of the input system.

V. CONCLUSION
We have presented a direct property estimation scheme, where a classical representation of quantum states is constructed with quantum reservoir processing and used for property estimation in the post-processing phase.Unlike existing techniques of shadow estimation that considers complex unitary ensembles and randomized measurements, our scheme explores the versatility of the quantum reservoir platform and requires only a single measurement setting.The pair-wise interacting reservoir nodes considered in our scheme results in minimal quantum hardware for the DV reservoir, and minimal training cost that depends only on the input system.The sample complexity has been rarely addressed in previous works on quantum reservoir computing.In contrast, we have established a stringent performance guarantee regard-ing the number of samples to be processed by the reservoir network.Furthermore, our scheme can be naturally extended to higher-dimensional systems and hybrid systems with non-identical local dimensions.We complement the theoretical results with diverse applications.For future research, reservoir computing is suitable for temporal pattern recognition, classification, and generation [37,78,79].While we have proposed a method for analyzing the statistical fluctuation of quantum reservoir outputs in property estimation, it is important to do so in temporal information processing tasks such as temporal quantum tomography [80] and nonlinear temporal machine learning [81].Also, it is important to further devise an operational toolbox for optimizing the reservoir networks regarding the sample complexity and investigate whether the efficiency lower bound for local measurements [5,11] can be achieved on the quantum reservoir platform.Moreover, there has been extensive research conducted on the learning abilities of quantum reservoir networks in tasks ranging from quantum to real-world problems [82][83][84][85].These studies explore various network topologies and connectivities.While the pair-wise quantum reservoir construction is beneficial for property estimation, there lacks a general understanding of the role of topology and connectivity in quantum reservoir computing regarding sample complexity.Besides, quantum reservoir computing can be carried out with a single non-linear oscillator [86], so it is important to analyze whether such a quantum reservoir can perform shadow estimation as well.Another topic is the noise effect.There are inspiring discussions on the noisy training data [87], robustness in tomographic completeness [88] and benefits of quantum noise [89].In our scheme the influence of time independent system noise on T cancels out by a direct observation of Eq. ( 13).However, a complete discussion on the noise effect and its mitigation in quantum reservoir processing is important for future experiments.
Proof.The single-snapshot estimator Ω is a linear combination of the measurement results of mutually exclusive readout operators {ô i } with the corresponding weights {w i }, i.e., Thus, the expectation of Ω is the inner product between W and X.The variance of Ω is where ⊙ represents the Hadamard product.Thus, where we have neglected the term −(W X) 2 .The inequality saturates when Tr(Oσ) = 0.
Since we have W = ⟨⟨O|T −1 and X = T |σ⟩⟩, we can directly compute the factor F res .Define an operator B such that then where || • || ∞ represents the spectral norm.Also, it is worth noting that estimating identity operator doesn't affect the variance.Define Ω′ = (W + cW 1 1 ) X, where W 1 1 = [1, 1, . . ., 1] is the weight operator for the identity observable, and and by construction i xi = 1, we have Var( Ω′ ) = Var( Ω).We observe that Ω0 results in a tighter worst-case variance upper bound on average in fidelity estimation of random pure target states.Also, for properties that only differs by a product factor c, With these properties, we note that if a reservoir setting works relatively better for overlap estimation of a target state σ, one could expect that it works better for estimating arbitrary linear properties that can be decomposed as An arbitrary normal matrix is related to a density matrix by Eq. (B9).Thus, a well trained reservoir setting for overlap estimation of random target states is applicable to the future estimation of multiple unknown normal matrices.In order to analyze the search of good reservoir parameters, we fix the reservoir parameters {J, P i , E i } and check the dependency of reservoir performance on evolution time.For the evolution time from t = 0 to t = 5.5, the variance upper bound of single-qubit state fidelity estimation averaged for 300 random target pure states varies, as shown in Fig. 8.It can be observed that the average variance upper bound is stable around some local minimums in a short time window.
Furthermore, we could analyze the efficiency scaling of estimating k-local tensor product observables with the help of pair-wise reservoir dynamics.Here we present the detailed proof for Theorem 2. W ⊙ W X is the worst-case variance upper bound for the single-snapshot estimator.

Multi-snapshot estimator
The variance can be further suppressed by combining single-snapshot estimators with statistical methods.The median of means method is expected to reduce the effect of outliers [5,8], the efficiency of which is given by the following theorem.3. Process N sample i.i.d.copies of the unknown state σ with QRP.For each copy the reservoir evolution time is randomly chosen from {t k } with the probability distribution P PTM .Load the readouts {X i } to classical memory.
4. Calculate the estimated values { Õi } with the weights and readouts.
Numerical result shows that for random reservoir dynamics, the variance could be suppressed by PTM.For each target state, we use the training data collected at 2 different time points with a random reservoir initialization, and optimize the variance upper bound with PTM.After that, we could achieve a more efficient estimation performance, as shown in Fig. 9.
property of median estimator ensures that if each U-statistic estimator has a variance no larger than ϵ 2 /34, then Pr Ω(2) i − Tr(O i σ ⊗2 ) ≤ ϵ ≥ 1 − 2e −K/2 .(C26) Thus, for each Ω(2) i the confidence level is no less than 1 − δ/M .The union bound ensures that the over all confidence level for estimating M properties is no less than 1 − δ.From Lemma.5, we choose the size for each set as 272/ϵ 2 max i A (2) (O i , X t ).Consequently, the total number of copies consumed is N sample = 544/ϵ 2 ln(2M/δ) max i A (2) (O i , X t ).
FIG.1.Schematic illustration of predicting many properties with quantum reservoir processing.A source prepares an n-qubit quantum state σ which is taken as input by a pair-wise reservoir network with 2n nodes.For the i-th input copy, the local measurement operator on the j-th node is nj = (1 1 − ς z j )/2, and the readout is a bit string si ∈ {0, 1} n , corresponding to one of the N read readout operators.Equipped with the training data, an unbiased estimator Ω is constructed for the property O.

FIG. 2 .
FIG. 2. Comparison with the shadow estimation withPauli measurements[5].(top) Average variance upper bound of fidelity estimation.For each system size we randomly generate 1000 pure states as the target states.The ratio of the average variance upper bounds for QRPE and shadow estimation remains almost constant as the number of qubits grows.(bottom) The average number of measurement settings invoked to reach a given precision with confidence level 0.90.For shadow estimation, we assume the variance in single-snapshot estimator equals a tenth of the shadow norm on the top figure.While QRPE requires only a single setting, the standard shadow protocol with Pauli measurements[5,49] requires exponentially increasing settings to reach an arbitrary precision.The reservoir setting we choose for qubit system is given by J = −0.41,P1 = 4.0, P2 = 1.3, E1 = 0.71, E2 = 0.46 and t = 1.The unit for the Plank's constant is meV•ps.

FIG. 3 .
FIG. 3. Observation of entanglement sudden death.(top): estimation of WGME.(bottom): estimation of WME.We use 6000 measurements for each {q, κt}, and the maximal estimation error for estimating GME and ME witnesses are 0.11 and 0.13 respectively.It can be observed that under a small dephasing noise, both the GME and ME witnesses vanish at a finite time for some initially entangled states, which indicates the emergence of ESD.

FIG. 4 .
FIG. 4. Numerical experiment for estimating Pauli observables of randomly chosen one, two, three and four-qubit input states.Each data point represents the average error of all the k-local Pauli observables of a random k-qubit input state.The solid marks represents the average performance of 30 independent experiments, while the hollow marks represents a single experiment.The subfigure shows the variance upper bound Fres, averaged over the k-local Pauli observables.Numerical results indicate that the worst-case sample complexity of estimating Pauli observables increases exponentially with the locality.

FIG. 6 .
FIG. 6. Saturation of the second Rényi entropy S2(ρA) of a small region A consisting of NA qubits.The solid curves represent S2(ρA) of a fourteen-qubit state generated by a VQE circuit with depth p, which approaches the page entropy represented by the dashed line.For reservoir estimation at each circuit depth we use 2000 independent snapshots.The shaded region is the fluctuation of estimated value in 10 independent estimations.

FIG. 7 .
FIG. 7. The virtual distillation of maximally entangled states mixed with depolarizing noise in 2 × 3 and 3 × 3 dimensional systems.The dashed curves represent the fidelity of the distilled states with m = 2, while the solid curves represent the fidelity of physical states with m = 1.The shaded region is the range of fluctuation of the estimated value in 10 independent estimations.For each reservoir estimation we use 10000 independent measurement readouts.The reservoir setting we choose for qutrit systems is J = 0.9, P1 = 2.1, P2 = 1.1, E1 = 1.1, E2 = 0.4, α1 = 0.6, α2 = 0.7.

Lemma 3 (FIG. 8 .W 1 1 ,
FIG.8.Average variance upper bound of single-qubit state fidelity estimation that varies with evolution time from t = 0 to t = 5.5.The parameters J, Pi, Ei are the same with main text.At each evolution time, the variance upper bound for pure state fidelity estimation Fres is averaged over 1000 random target states.While there are multiple local minimums with different slopes at neighboring points, the performance of reservoir parameters remains stable for a short time period at some local minimums.

FIG. 9 .
FIG. 9. Variance reduction via PTM.For each random test we use the training data collected at t = 1 and t = 10 of a two-node reservoir with parameters J, Pi, Ei uniformly generated at random from [0, 5].The variance upper bound for overlap estimation Fres is averaged over 300 random target states.The green bars uncovered by blue bars indicate the improvement in average variance upper bound achieved by PTM.
Each training step starts with inputting a training state |φ k ⟩⟨φ k | into the reservoir to reach the initial training state ρ(0) = |φ k ⟩⟨φ k | ⊗ [|0⟩⟨0|] rest , which then evolves to ρ(t) at time t.From sufficiently many measurement results of each input training state, we estimate the expectation value ⟨ô i (k)⟩ of the readout operator.The training data is stored by arranging ⟨ô i (k)⟩ into a column vector X|φ k ⟩⟨φ k | .In this way, we collect readout vectors X|φ k ⟩⟨φ k | corresponding to all the training states |φ k ⟩ for k = 1, 2, . . ., N train .
Here we consider the training data to be accurate, and present results that account for statistical noise occurring outside of the training phase.The reservoir dynamics are initialized by setting the parameters J, P 1,2 , E 1,2 and the evolution time t.
The optimization over probability distributions aims to reduce the average variance upper bound of single-qubit state overlap estimation, see Appendix B for more details.The QRPE protocol for linear functions is as follows:1.Perform a one-time estimation of the training states of a two-node reservoir and load the training data {X p } to classical memory.2.Given properties{O i }, calculate weights {W i } with the training data.Obtain the worst-case variance upper bound max i ||B i || ∞ .Calculate N sample with the given confidence 1 − δ and additive error ϵ. 3. Process i.i.d.copies of the unknown state σ with the quantum reservoir network.Load N sample snapshots {X i } to classical memory.
, we could utilize random settings by probabilistic time multiplexing (PTM), where the pair-wise training data X p (t k ) are collected at N time different time points {t k | k = 1, 2, . . ., N time }, and the single-snapshot estimator is constructed by measuring reservoir nodes at time t k with probability p k , where k p k = 1. 4. Calculate the estimated values { Õi } with the weights and snapshots.