Online quantum time series processing with random oscillator networks

Reservoir computing is a powerful machine learning paradigm for online time series processing. It has reached state-of-the-art performance in tasks such as chaotic time series prediction and continuous speech recognition thanks to its unique combination of high computational power and low training cost which sets it aside from alternatives such as traditionally trained recurrent neural networks, and furthermore is amenable to implementations in dedicated hardware, potentially leading to extremely compact and efficient reservoir computers. Recently the use of random quantum systems has been proposed, leveraging the complexity of quantum dynamics for classical time series processing. Extracting the output from a quantum system without disturbing its state too much is problematic however, and can be expected to become a bottleneck in such approaches. Here we propose a reservoir computing inspired approach to online processing of time series consisting of quantum information, sidestepping the measurement problem. We illustrate its power by generalizing two paradigmatic benchmark tasks from classical reservoir computing to quantum information and introducing a task without a classical analogue where a random system is trained to both create and distribute entanglement between systems that never directly interact. Finally, we discuss partial generalizations where only the input or only the output time series is quantum.


Introduction
Tasks where one time series need to be transformed into another include time series forecasting 1,2 , pattern generation 3,4 and pattern recognition [5][6][7][8] .In online time series processing both the given data and desired transformed data are functions of time, which separates it from approaches such as first recording the data and later processing it.Instead, the objective is to realize the time dependent function which for a given timestep and input time series up to that step returns the corresponding element of the output time series.Such tasks are also known as temporal tasks.When successful, online time series processing facilitates, e.g., the processing of arbitrarily long sequences of data since the inputs are continuously processed into outputs.This is possible in particular for tasks that can be solved by so called fading memory functions, which are functions well approximated by continuous functions of only a finite number of past inputs 9 .Under typically mild conditions the input time series can be used to drive random dynamical systems such that their internal variables become such fading memory functions, which can then be combined to approximate the desired ouput by training a simple, even linear readout function.This is known as reservoir computing (RC), which is a powerful approach to solving temporal tasks thanks to a remarkably low training cost 10,11 combined with state-of-the-art performance 12 .Furthermore, classical or quantum physical systems are also amenable to be used as the dynamical system [13][14][15] , paving the way to harvesting computational power from essentially random physical systems with fading memory and complex dynamics.In RC such systems are usually called reservoirs.
After recent seminal works investigating the suitability of the transverse-field Ising model for RC purposes 16,17 , there has been a surge of interest in the quantum case in particular.Indeed, the initial model has been refined and analyzed in several ways [18][19][20][21] , whereas new proposals have introduced RC based on quantum circuits 22,23 , NMR systems 24 and continuous variable quantum systems 25 .The results have been promising, suggesting that both in the discrete and continuous variable case quantum reservoirs may have an advantage over their classical counterparts in terms of how rapidly the potential reservoir performance improves with size 16,20,25 .One of the biggest hurdles is in fact the extraction of the classical output from the quantum systems, since not only does a single measurement reveal only a tiny amount of information about a quantum system in an unknown state, it also alters the state and therefore competes with the inputs in driving reservoir dynamics.For certain special systems, such as NMR systems, an enormous amount of copies of the reservoir are naturally available which has been proposed to allow one to bypass the measurement back-action problem when collective input injections and measurements can be carried out 24 .In general, repeatedly initializing and subsequently measuring a quantum system extracts classical information out of it-for example, a value of an observable-however such an approach is far from ideal for time series processing for two reasons.Firstly, carrying out the repetitions anew for every element in the output time series one wishes to learn introduces severe overhead, and secondly, it is hard to imagine how such a protocol can run in an online mode, continuously producing elements of the output time series.All in all, output extraction is a major challenge in exploiting the potential of quantum reservoirs for classical time series processing.
Here we lay down an alternative, RC inspired path that can fully harness the quantumness of the reservoir while largely sidestepping the measurement problem.Namely, we introduce online time series processing with random fading memory quantum systems where both the input and desired output time series consist of quantum information.Such temporal quantum tasks are by construction impossible for a classical reservoir, whereas no measurements are required after the reservoir has been trained since the output can remain quantum.Specifically, we consider random networks of interacting quantum harmonic oscillators as the reservoir.Taking inspiration from RC, we train only the interaction Hamiltonian between the network and the carriers of input information.The main difference with RC is how the output is formed; in the former it is a trained function of reservoir observables, here the output is imprinted directly on the quantum systems acting as carriers of data.This affects also the training process as will be seen.
We illustrate the possibilities of the proposed model with three different temporal tasks.The short term quantum memory (STQM) task is the quantum analog of the short term memory task 26 commonly used as a benchmark task in classical RC-the objective is to recall past inputs that are no longer available using the memory of the reservoir.Another common task is channel equalization task 27 , where the input time series is transmitted through a noisy nonlinear channel that also mixes the time series with various echoes of itself, and the objective is to recover the original time series from the distorted one.Here we generalize it to inverting the transformation caused by a quantum channel.Despite being generalizations of a classical task it will be seen that the quantum cases have notable differences.Furthermore, we introduce a task without a classical counterpart which we call the entangler.Here the objective is to create entanglement between different initially uncorrelated systems by letting each of them in turn interact once with the reservoir but never between each other.We find that all these tasks are possible to solve using random untrained networks of interacting quantum harmonic oscillators; remarkably, not even the network initial state needs to be controlled.Finally, we briefly discuss partial generalizations, i.e. cases where only the input or only the output time series is quantum.

The model
In RC a reservoir is a dynamical system that can be steered by an input time series to a trajectory in its state space determined by the inputs alone, i.e. its internal variables become completely determined by the input history at the limit of many inputs.If the variables can be monitored then the response of the reservoir at different timesteps can be post-processed to achieve a desired transformation from the input time series to an output time series.Importantly, for sufficiently complex reservoirs nontrivial transformations can be achieved by cheap post-processing, such as a linear combination of the variables.
Here the reservoir is a network of N unit mass quantum harmonic oscillators interacting with springlike couplings.Such units are used that h = 1 and k B = 1.Let p = {p 1 , p 2 , . . ., p N } and q = {q 1 , q 2 , . . ., q N } be the vectors of momentum and position operators of the oscillators.The reservoir Hamiltonian H R is where the diagonal matrix ∆ ω holds the oscillator frequencies ω = {ω 1 , ω 2 , . . ., ω N } and the symmetric matrix L has elements Here g i j ≥ 0 are interaction strengths between the reservoir oscillators.Aside from oscillator frequencies there is a one-to-one correspondence between H R and weighted simple graphs.Indeed, L can be interpreted as the Laplace matrix of such a graph, and a given graph with a Laplace matrix L defines H R through Eq. (1).We consider temporal quantum tasks (analogous to temporal tasks in RC), which we define in this work as follows.The input time series s = {. . ., ρ I m−1 , ρ I m , ρ I m+1 , . ..} consists of quantum states ρ I m where m indicates the timestep.In general these can be states of multimode continuous variable quantum systems, however we assume that apart from their states the systems are identical, i.e. each system has the same Hamiltonian H S .The input time series is processed by the reservoir into output . .} by letting each system in turn interact with the reservoir for some time ∆t according to an interaction Hamiltonian H I coupling every reservoir oscillator to every subsystem.The order of interactions is given by the timesteps.Consequently o is the image of s and reservoir initial conditions under a transformation induced by the full Hamiltonian H = H R + H S + H I and the interaction time ∆t.In a temporal quantum task we attempt to realize a given transformation from s to o in this way.Besides the uncorrelated case one may also consider correlations between the systems at different timesteps, and we will return to this point later.
In the special case where H S consist of M unit mass quantum harmonic oscillators and the interactions in H I are springlike couplings, H has the same general form as H R .The transformation induced by H and ∆t on the operators of the reservoir and input system is now linear and can be given in terms of a symplectic matrix S. Let x R k be the form of the reservoir operators after k-th input has been processed, let x S k be the operators of the k-th input and let x O k be the operators of the k-th output.Now where the symplectic matrix has been divided into blocks such that A is 2N × 2N and D is 2M × 2M.By iterating this equation we immediately get the form of both the reservoir and input modes for some timestep m: where x R 0 is the initial form of the reservoir modes.The form of the output modes for some timestep m as a function of x R 0 and input history is then where the first line is exact and the second line an approximation which holds when the spectral radius ρ(A)-not to be confused with quantum states-is less than 1 and enough inputs have been processed.Given that the equations of motion are conveniently expressed in terms of the operators, in the rest of this manuscript we will simply write s = {. . ., x O m+1 , . ..}.When ρ(A) < 1 the output time series o becomes independent of the initial conditions at the limit of infinitely long input history, which is known as the echo state property 1 in the RC literature.In fact, it can be shown 25 that satisfying the spectral radius condition gives the reservoir also the so-called fading memory property 9 , which guarantees that x R m and therefore x O m become well-approximated by a continuous function of only a finite number of past inputs at the limit m 1.This not only ensures that the initial conditions can be ignored but also prevents any physical quantities from diverging: the reservoir state never leaves the state space as long as all input states are physical, not even in the limit m → ∞.At variance, if ρ(A) ≥ 1 then, e.g., reservoir excitations may diverge.Finally, in the special case where A is nilpotent for some index n there is sudden death of reservoir memory where x R m is a function of exactly n previous inputs.As the index of a nilpotent matrix is always at most its order 28 , n ≤ 2N.

RC inspired quantum time series processing
The system given in Eqs.(3) can be harnessed for RC by considering only x R m when it is assumed that the states in s are in fact functions of the elements of a classical time series 25 .When ρ(A) < 1 there is fading memory and the observables of the reservoir become well-approximated by continuous functions of only a finite number of past inputs.Different transformations to a classical output time series can then be accomplished by training a simple function of the reservoir observables such as first moments, second moments or covariances of x R m .The resulting RC can be analyzed with contemporary RC theory since the latter is agnostic to the mechanism that creates the functions of the input.Importantly, H R is not trained and can be random since ρ(A) < 1 can typically be achieved by just tuning ∆t.
The scheme for using the reservoir to process temporal quantum information instead is shown in Fig. 1a.Here s itself will be the input while o consists of quantum information encoded in x O m .Thanks to fading memory, each x O m is completely determined by x I i where i ≤ m and the initial reservoir state can be ignored.To achieve different transformations s → o, the matrix S induced by H = H R + H S + H I and ∆t must be changed while preserving ρ(A) < 1.Whereas the number of parameters in H is proportional to (N + M) 2 , we take an approach inspired by RC and attempt to achieve online quantum time series processing by only training the NM interaction terms in H I , leaving both H R and H S fixed.We will provide strong numerical evidence that this is in practice enough to succeed in many different temporal quantum tasks, even with random H R .
The similarity of this result with the fact that in RC with fading memory systems it is enough to only train the final weights in the readout layer is uncanny, since the situations are quite different.Indeed, unlike in RC here reservoir dynamics cannot be separated from training since tuning H I changes the dynamics.We are in fact not aware of any theoretical results that could be used to explain the phenomenon.In the following we provide a brief overview of the general purpose training process.See Methods for further details.First the input time series s is divided into three phases: the preparation, training and test phases.The role of the preparation phase is to get rid of the influence of the reservoir initial state.During the training phase the performance is checked using a cost or objective function, which varies depending on the task but for a fixed input and reservoir is completely determined by H I .The specific forms will be introduced along with the respective tasks.The interaction Hamiltonian is varied to optimize the function using a simple stochastic function optimizer.Spectral radius condition is enforced by providing the optimizer as initial points only such H I that the condition is satisfied; during minimization points that violate the condition can be expected to perform worse and are therefore discarded.Unless the time between inputs ∆t is fixed by the task, training may be repeated for different choices of ∆t.In the test phase trained H I and best ∆t are used and reservoir output is collected to check the performance using a task dependent figure of merit.In this way the trained reservoir is exposed to new input and must be able to generalize beyond the specific inputs in the training phase to succeed.These are the results shown in the figures.

Parameter values used in numerical experiments
In all numerical experiments reported throughout the section we consider as the reservoir random completely connected networks of N identical oscillators with a bare frequency ω 0 = 0.25.Each coupling strength between the reservoir oscillators is g i j ∈ [0, 0.2], chosen uniformly at random.The M input modes are are also such oscillators.Although in principle having ρ(A) < 1 is enough for fading memory, in practice we impose the limit ρ(A) < 0.99 to avoid issues with finite precision numerics.The influence of the initial state of the reservoir is washed out during the preparation phase and is therefore irrelevant, however in practice we use the ground state of H R .For each different case we show results of 100 random realizations of the reservoir and when applicable, other quantities such as the input time series s.In all cases the lengths of preparation, training and test phases are 40, 80 and 40, respectively.The values of the time between inputs ∆t and the input states themselves vary and will be reported along the tasks.

Short term quantum memory task
The short term memory task is a paradigmatic task in classical RC where the input s k at some timestep k is a real scalar or vector, and the target is s k−τ where τ is the delay.Checking the performance in this task for different delays is often used to gauge how much linear memory a reservoir has.The short term quantum memory (STQM) task is its direct generalization where s k is a state of a quantum system.Specifically, The objective is to achieve the transformation s → o defined by Eq. ( 5) by training ∆t and Hamiltonian H I .H R is assumed to be random but fixed and x R 0 can be arbitrary.Unlike in the classical case where the amount of information grows linearly with input size, here the growth is more rapid as the reservoir must be able to delay also the correlations and entanglement between different input modes in x I k .We mention in passing that the process tomography of such a delay map in the discrete variable case was considered in Ref. [29].
It should be pointed out that for single-mode input states the task can in principle be done exactly for any delay τ by concatenating N = τ single oscillator reservoirs, which can be compared to deep RC.For τ = 1 the states of the reservoir and input must be swapped; the required interaction strength and time can be solved analytically.Clearly using a single reservoir oscillator with double the interaction time solves the task for τ = 0 since the states are swapped twice, while using a sequence of swaps with N different single mode reservoirs achieves a delay τ = N.For M > 1 and τ > 0 the task can be done by using N = Mτ non-interacting reservoir oscillators, provided that the input systems are likewise non-interacting.Here we show that the task can in fact be solved by using a single random and fixed reservoir and letting the input interact only once with it, which can be compared to ordinary (shallow) RC.
The simplicity of this task makes it amenable to a special purpose training procedure.Indeed, from Eq. ( 4) the following conditions to solve the task can be observed: When satisfied the contributions from the incorrect timesteps are suppressed while the contribution from the correct one is enhanced.The conditions in Eqs.6 for the full symplectic matrix induced by all three Hamiltonians must be satisfied by training only one of them, H I .In practice the cost function given by Eqs.11 in Methods is minimized to train the reservoir.Additionally, training of H I is repeated for ∆t = 2π/ω 0 , 4π/ω 0 , . . ., 8π/ω 0 and the best value is chosen for testing phase.Notably, training is input state independent; although the training phase in s could therefore be omitted, it is kept for the sake of consistency.It should be stressed that since the task is linear in the involved modes for any state, if training is successful the reservoir can delay also non-Gaussian M-mode states.
Results of numerical experiments are shown in Fig. 2. The input time series s consists of random zero mean M-mode Gaussian states generated by acting with a random symplectic matrix on a thermal state of non-interacting oscillators.For details and parameter values used see Methods; it should be stressed that in general this creates correlated states.The reservoir parameters are as specified previously.Different values of N, M and delay τ are considered, whereas the figure of merit is the (Uhlmann) fidelity averaged over the test phase, possible to calculate in closed form for arbitrary Gaussian states applying, e.g., results of Ref. [30].
In panel a reservoir size N and input size M have been fixed and the delay τ is varied.The performance is excellent especially for small delays τ = 0, 1, 2 where the median values are F > 0.999999, F > 0.9999 and F > 0.964, respectively.Even the worst possible performance is decent for small delays, and it is conceivable that if in their case training was repeated with a different seed for the random number generator the performance might be improved significantly.In panel b the delay is fixed and N and M are varied; the diagonal part of the array shows the case N = M.As expected, the reservoir struggles significantly with the task when N < M whereas for N > M the performance is mostly good.Remarkably, a random reservoir can be trained to delay also multimode states by only tuning H I and ∆t.

Quantum channel equalization
This task is inspired by the channel equalization task where the input time series is distorted by a transmission through a classical channel with fading memory and the objective is to invert the transformation by the channel and restore the original time series.The quantum counterpart is presented in Fig. 1c, where the original time series consists of states of quantum systems with operators x I that are transformed by an interaction with some fixed but random system with operators x Ch , which constitutes the channel.We assume the channel to induce a linear transformation of the modes given by some symplectic matrix, i.e. the relevant Hamiltonians are quadratic, and furthermore that it has fading memory, which implies that x Ch is   well-approximated by a function of a finite number of past inputs.The transformed operators of the input system are denoted by x D .A reservoir with operators x R is trained to recover the original time series from x D .The equations of motion can be recast in the same general form as before as explained in Methods.
It should be stressed that simply training the reservoir to perform the inverse of the channel symplectic matrix will not work since the channel acts on x Ch and x I but the reservoir symplectic matrix S acts on x R and x D .The transformed input x D k for some timestep k in general does not contain full information of any of the previous inputs because part of the information is in x Ch and in general also in the correlations between the channel and the distorted input.Since it is well known that unknown quantum states can neither be cloned nor amplified the task as given is in fact impossible-in stark contrast with its classical counterpart.To have any hope of success, the task must be modified.
Here we will do this using techniques inspired by classical RC.Instead of a single copy of the input state we will transmit a product state of s copies, but still require the reservoir to distill only a single copy of the original input, thus giving the reservoir additional quantum information to work with.We call this spatial multiplexing at order s.Additionally, we transmit the same product state m times, one copy after another, requiring the reservoir to output the original input only after the m-th copy.We call this temporal multiplexing at order m.More formally, Here the timesteps are to be understood as the points where we switch from one set of identical copies to another, as determined by the original unaltered input time series s.As a final remark before moving on, if s is given but unknown neither spatial nor temporal multiplexing can be used, making the task again impossible.It must be assumed that there is a source that directly generates states according to some s and some m, or alternatively that the states are known to the sender who may then prepare the copies.It may be asked if some other modifications could help solve the task even for unknown s, but this is outside the scope of present work.With the analysis complete for now, we move on to numerical experiments.The input consists of single mode zero mean Gaussian states.The channel has C = 2 oscillators and the Hamiltonian has the same general form as the reservoir Hamiltonian.This system is taken to interact with the inputs such that the spectral radius of the relevant block in the symplectic matrix is at most 0.95.For the reservoir N = 3.The input time series consists of random zero mean Gaussian states like before in the STQM task, however for simplicity we focus only on the single mode case where M = 1.Unlike before, the time between   inputs is taken to be fixed by the channel and is therefore kept at a constant value which was chosen to be ∆t = 1.5π/ω 0 .Results are shown in Fig. 3.
In panel a spatial and temporal multiplexing are considered separately.At s = m = 1 both reduce to the original formulation of the task, which was already concluded to be unsolvable.Still, performance exceeds that of random guessing.Performance increases quickly with s and slowly with m.Indeed, already at s = 2 results tend to be better than at m = 5.That being said, unlike increasing m increasing s increases the number of terms in H I which, e.g., makes training slower.In panel b spatial and temporal multiplexing are considered together.Curiously, performance does not always increase with m and s.For example, both m = 3, s = 4 and m = 4, s = 3 lead to a better performance than m = 4, s = 4.Moreover, spatial multiplexing alone achieves a performance not too far off from the best case.

Entangler
In this task the inputs s k are taken to be uncorrelated single mode systems in the vacuum state whereas the target time series consists of entangled states.Although more complicated patterns can be envisioned, here we focus on entangling systems with a fixed delay τ.That is to say the goal is to entangle s k with s k−τ by training H I .For τ = 1 the target is then a chain of systems with nearest neighbor connections where the connections are entanglement, for τ = 2 a chain with next nearest neighbors connected and so on.Importantly, we assume that only one input interacts with the reservoir at any given timestep and all the others are unavailable, and furthermore assume that there are never any direct interactions between the systems in s.
In fact, if we imagine that the systems are periodically generated by a source of vacuum states then the input s k+1 does not even exist at some timestep k.A system or a device that can solve the task can turn the source of uncorrelated states into one of entangled states.
Much like in the STQM task, in certain special cases and allowing for time-dependent H I the task can be achieved analytically and exactly, as depicted in Fig. 1d where only the case τ = 1 is considered for simplicity.At each timestep the reservoir and the ancilla are first entangled using an interaction Hamiltonian of one form and then their states are swapped using a different interaction Hamiltonian; in both cases one can analytically solve what the precise form of the interactions and interaction times must be.In fact, since the operations commute the order does not matter.Every application of the entangling gate creates a link in Fig. 1d, which is later swapped to the next input system.The role of the reservoir oscillator is to provide short term quantum memory.Without it, the task becomes impossible.
Here we solve the task with the RC inspired approach by training the time-independent H I to maximize the entanglement, as quantified by logarithmic negativity 31 between s k and s k−τ .Like in the STQM task, training is repeated for ∆t = 2π/ω 0 , 4π/ω 0 , . . ., 8π/ω 0 and the best value is chosen for testing phase.To succeed the reservoir must simultaneously create entanglement and re-distribute the quantum information correctly as explained previously.Results are shown in Fig. 4  panel a the reservoir size is fixed to N = 3 and the delay τ is varied, and the logarithmic negativity averaged over the systems in the test phase is shown.The logarithmic negativity achieved for the shortest delay is between 0.35 and 0.4, which corresponds to that of a twin beam state with two-mode squeezing parameter 0.175 ≤ s ≤ 0.2.Performance decreases slowly up to N = τ, but then collapses for delays τ > N, bearing a striking similarity with panel a of Fig. 2. One may interpret this as the reservoir being able to remember up to N single mode states before running out of memory.The same behaviour can be observed also in panel b where both N and delay τ are varied and median performances are shown.

Partial generalizations
In all previously introduced tasks both the input and the output time series consist of quantum information, however one may consider partial generalizations where one of them is still classical.Here we briefly illustrate the possibilities with two simple examples.
If the output is quantum but s is classical, say, a time series of systems in thermal states, one may follow the framework used previously.As an example task we consider predictive quantum state preparation where the reservoir is trained to prepare a given quantum state-here, squeezed vacuum-based on future classical inputs.For arbitrary s this is of course impossible, but if s is at least approximately predictable then the task can be in principle solved.Here we consider the Santa Fe chaotic time series, a dataset recorded from a far-infrared laser in a chaotic state 32,33 often used to benchmark the predictive power of classical reservoirs.Specifically, we normalize the Santa Fe time series and consider as s single mode thermal states such that the number of thermal excitations (n th ) k follows the normalized time series and the target is a squeezed vacuum state with a squeezing parameter r k = (n th ) k+a where a ≥ 0 is the advance, or the number of timesteps in the future the reservoir must be able to predict.
Results are shown in Fig. 5a, where the average fidelity between the target state and the actual output state is shown for different values of the advance.There is very little spread for most values since only the reservoir is randomized between different realizations.Interestingly, there is an abrupt change in behavior when the advance a exceeds the number of reservoir oscillators N.Even then the fidelity remains decent but there is considerably more spread in the performance.
In the opposite case where output is classical information, say, about the properties of the states carried by the input systems in s, the approach where the classical output is formed from reservoir observables can be used, as outlined previously in the Model section.Let σ (x R k ) be the reservoir covariance matrix at some timestep k.It can be shown 25 that which holds for any number M of input modes.In principle, the elements of σ (x R k ) can be estimated by performing measurements on multiple copies of the reservoir that has processed identical inputs s, which however introduces substantial overhead.Analyzing just how much overhead is incurred is beyond the scope of this work; in what follows, it is assumed that the exact values of the elements of σ (x R k ) are available.That being said, since the target is some function of reservoir observables the full Hamiltonian can remain constant, decoupling training from the dynamics.Indeed, once the elements are available multiple trained functions can be used to estimate a number of different features of s.
As an example task, we consider as input random single mode Gaussian states and as target det(σ (x I k−τ )), or the determinant of the single mode covariance matrix for some delay τ ≥ 0. This is an important quantity that, e.g., completely determines the purity, amount of thermal excitations and the von Neumann entropy of the state.Below we given an overview of the conditions under which this task was simulated; for full details, see Methods.
We consider a reservoir of size N = 20 and consider M = 10 input modes such that the input is in a random product state of identical single mode states.Furthermore, the Hamiltonian is such that only two reservoir oscillators interact with a single input mode, and there are no interactions between these triplets of two reservoir oscillators and a single input oscillator.The interaction strengths are random but fixed and the spectral radius condition is satisfied by tuning ∆t.One may observe from Eq. ( 8) that σ (x R k ) is linear in σ (x I k−τ ) for any delay τ unlike the determinant, however this problem can be overcome by considering trained linear combinations of products of pairs of elements of σ (x R k ).Here the output is a trained linear combination of products of distinct pairs of the first row of σ (x R k ), with training carried out as in Ref. [25].Finally, unlike elsewhere, we consider preparation, training and test phases of length 500, 2000 and 500, respectively.
Results are shown in Fig. 5b, where the the normalized mean squared error (NMSE) between the actual von Neumann entropy to that computed from the determinants estimated by the reservoir is shown.As can be seen, the NMSE is very small for all considered delays, suggesting an excellent agreement between the actual and predicted value.

Discussion
In this work we have introduced a RC inspired model for online processing of time series consisting of quantum information.Importantly, we have found that just with a judicious choice of the interaction Hamiltonian random instances of the model starting from any initial state can solve a variety of different tasks with high performance.We have also briefly illustrated the possibilities of partial generalizations of classical temporal tasks to cases where either the input or the target time series remains classical.Looking at the bigger picture, it is interesting to compare and contrast two distinct situations: when the output time series is to be quantum, and when it is to be classical.
In the former case the output extraction problem hindering previous related work vanishes but engineering freedom is preserved: control of only a small subset of all parameters is sufficient for high performance.Furthermore, if the reservoir Hamiltonian is random but known then measurements are not needed even in the training stage provided one can simulate the dynamics.For the considered model in particular an unknown Hamiltonian can be probed first 34 .It can be imagined that a classical RC augmented with a state preparation mechanism could emulate the case where input is classical, but otherwise there is genuine quantumness: in the single shot case it is clear that no classical RC can emulate its quantum counterpart since the input data would first have to be transformed to classical information.That being said, training the interaction Hamiltonian is in general somewhat costly even when the dynamics can be simulated.If simulation is not possible the cost or objective function must be estimated with measurements, which should be expected to be a very challenging optimization problem in its own right by comparing with, e.g., variational quantum algorithms [35][36][37][38][39] .Finally, another advantage of RC is lost in multitasking where the same reservoir processing the same input can simultaneously solve many different tasks by using differently trained readout functions.In the quantum case any attempts to multitask will inevitably affect performance because there is only so much quantum information for forming the output.
When instead the output is classical, multitasking is possible and the training cost is minimal.The challenges are two-fold: the output extraction problem and pinpointing what role exactly quantumness plays aside from providing a larger state space.Moreover, even if the output extraction problem can be solved, the specific way it is solved may dictate what quantum systems are ultimately suitable.As we demonstrate here with von Neumann entropy detection, the case where input is quantum might be of particular interest however, since thanks to the memory and multitasking a plethora of information concerning multiple past input states can be distilled even if the reservoir observables are known only for one or few timesteps.This may be compared to recent proposals where information of only a single quantum state is extracted with the help of a larger quantum system and supervised machine learning [40][41][42] .
Indeed, the model proposed here and the results have created a fertile ground for further work in the direction where at least one of the time series is quantum.To the best of our knowledge there is currently little work on such temporal tasks, however the inverse problem of performing tomography of an unknown temporal quantum map has been recently considered in spin system 29 .Comparisons may also be made with a recent proposal to train quantum system to induce quantum gates between qubits 43 ; its temporal generalization might consider gates between inputs at different timesteps, for example.Indeed, temporal quantum tasks could be tackled also in the discrete variable case.There is also room for further improvements for the introduced model by considering for example how training also the time between inputs can affect the performance or considering the case of non-Gaussian states or operations, which can be expected to lead to nonlinear memory 44 where x O m can be nonlinear in x I i for i ≤ m.One may also consider the prospects of a proof-of-principle experimental implementation since the general form of the reservoir Hamiltonian can in principle be realized in a multimode optics platform 45 .

Generation of random zero mean Gaussian states
In the single mode case where M = 1 the states may be parameterized in terms of the thermal excitations n th , magnitude of squeezing r and phase of squeezing ϕ.For displacement we consistently use α = 0, leading to the input first moments to vanish-the state is now completely characterized by its covariance matrix, which reads which is a covariance matrix of a squeezed thermal state.
In the multimode case where M > 1 the state is essentially parameterized by M thermal excitations, each independently and uniformly distributed, and M squeezing parameters, also independently and uniformly distributed, in the following way.We begin from the product state of M single mode thermal states, each with their own thermal excitations.Then we act with a random basis change, apply single mode squeezing of the position to all the modes with random magnitudes, and finally act on the resulting state with another random basis change.The random basis changes are built from Haar random M × M unitary matrices; let such a matrix be U. Then by construction is orthogonal and also a symplectic matrix w.r.t. the chosen ordering of operators.
For both STQM and channel equalization tasks we have chosen the intervals to be n th ∈ [0, 10] and r ∈ [0, 1].In the single mode case ϕ ∈ [0, 2π].These intervals also apply to the input states used for von Neumann entropy detection task shown in Fig. 5b.In the entangler task the input states are always single mode vacuum states.

Cost function minimization
A simple stochastic function optimizer called differential evolution (DE) is used.It treats the cost function as a black box, allowing it to, e.g., attempt to optimize functions where gradients (1st derivatives) or hessians (2nd derivatives) either do not exist or are not practical to calculate.Specifically, the implementation of Wolfram Mathematica 11.2 is used, which is described in Ref. [46].Here we give an overview of the method and the parameter values used; for full details consult the reference.
DE iterates a population of points {x 1 , x 2 , . . ., x d }.At each iteration a new population is created from the old one as follows.For each x j in the old population, three other old points x w , x u and x v are chosen randomly and a point x s = x w + s(x u − x v ) is formed where s ∈ R is a parameter called scaling factor.Then a new point x new j is created by taking each element either from x j or x s with probabilities p and 1 − p, respectively, where the parameter p is called cross probability.Finally, the new point x new j replaces x j if f (x new j ) is better than f (x j ), where f is a given cost or objective function.The stopping criterion is met when both | f (x new j ) − f (x j )| and x new j − x j are sufficiently small.We initialize the population by generating 30NM points, each corresponding to different interaction Hamiltonian H I where each interaction strength g nm between some reservoir oscillator n ∈ {1, . . ., N} and some input mode m ∈ {1, . . ., M} is uniformly and independently distributed in g nm ∈ [0, 0.2] such that the spectral radius condition ρ(A) ≤ 0.99 is satisfied.In the event that some point does not satisfy ρ(A) ≤ 0.99 it is generated anew.
We consistently use a scaling factor of s = 0.05 and a cross probability p = 0.4, i.e. rather small shifts are used to create the shifted points x s and when forming x new j the elements are slightly more likely to be picked from x s .We settled for these values through a simple lattice search.All other settings use the default values listed in Ref. [46].

Cost function of the STQM task
The cost function is where • is the Frobenius norm and where with a slight abuse of notation we have indicated by H I ∞ the maximum coupling strength between a reservoir oscillator and the input oscillator(s).The point of the term 1/ H I ∞ is to prevent the training to converge to the trivial solution H I = 0.The factors 0.5 and 5 control the relative importance of minimizing the norm of D and achieving CA τ−1 B ≈ I; these values where chosen after some trial and error.While the function does not feature all of the relevant terms in Eqs.(6), numerical experiments suggest that including more terms leads to worse results.

Objective functions of the quantum channel equalization and entangler tasks
Unlike in the relatively simple STQM task, in these tasks there is no obvious way to derive conditions on the reservoir symplectic matrix.This is why the objective function is the task dependent figure of merit-fidelity between reservoir output and original input in channel equalization and the logarithmic negativity in entangler-during training phase.

Additional details about the quantum channel equalization task
Let us write down the transformations caused by the channel and the reservoir at some timestep k.The interaction between input and channel modes induces a symplectic matrix S .Its action on all of the relevant modes reads where S has already been divided into blocks such that A is C ×C and D is M × M. Nothing happens to the reservoir modes since there is no interaction between the reservoir and the channel.The reservoir processes where the intermediate form x D of the input modes has been eliminated.The dynamics now follows Eqs. ( 2) through (4) with the replacements that is to say the channel and the reservoir may be treated together as if they formed a new, larger reservoir.This simplifies the equations of motion and the simulation of the dynamics.Although one may now consider Eqs.(6) to solve the task, in practice the performance is very poor because only the reservoir blocks are controllable, hence the modifications of Eqs.(7).

Additional details about the von Neumann entropy detection task
Let ρ be a single mode Gaussian state.Then its von Neumann entropy is defined as S V (ρ) = −Tr(ρln(ρ)).It can be shown 47 that S V (ρ) = n th ln n th + 1 n th + ln(n th + 1) (16)   where n th is the amount of thermal excitations of the state ρ.This quantity in turn is connected to the determinant of the associated covariance matrix σ through Det(σ ) = (0.5 + n th ) 2 , ( which can be seen by direct calculation starting from, e.g., Eq. ( 9).In Fig. 5b the actual S V (ρ) is compared to that computed from the estimated determinant of the input covariance matrix using Eqs.( 16) and ( 17).

Figure 1 .
Figure 1. a A time series of quantum systems with operators x I is processed into another time series where the transformed operators are indicated by x O .Each system interacts one after another with a random network of oscillators-later called the reservoir-with operators x R .In general the reservoir state depends on the states of systems it has interacted with, in turn making x O a function of all previous x I .Different transformations can be achieved by tuning only the interaction terms, indicated by dashed black lines.b In the short term quantum memory task the objective is to transform x I at timestep m into x Iat timestep m − τ where τ ≥ 0 is a delay.This is possible when the sought state can be distilled out of the reservoir memory, stored in x R .c In the quantum channel equalization task the input is transformed by a random, uncontrollable system with operators x Ch into the distorted input x D .The reservoir is to recover the original input from the distorted one.d In the entangler task the reservoir is to entangle the systems in the time series.Like before the systems never directly interact and only one of them interacts with the reservoir at any given time.Here there is only one reservoir oscillator and only a part of the input and output time series is shown for simplicity, whereas entanglement is indicated by solid black lines.

Figure 2 .
Figure 2. Results for STQM task.Different reservoir sizes N, input sizes M and delays τ are considered, and for each different set of values 100 random reservoirs were used.a Delay τ is varied for fixed N and M and all results are shown with a box plot.The box plot shows the minimum (lower whisker), maximum (upper whisker), median (line between boxes) and the first and third quartiles (beginning of lower box and end of the upper box, respectively).The fidelity achieved by random guessing is indicated by the dashed horizontal line.b Delay τ is fixed while both N and M are varied.Here only the median value is shown.See text for details.

Figure 3 .
Figure 3. Results for the quantum channel equalization task.Reservoir, input and channel sizes are fixed to N = 3, M = 1 and C = 2, respectively, while the orders of spatial and temporal multiplexing of the input are varied.In the former product states of identical copies of the input are used instead of single mode states, and in the latter, identical copies of the input are injected sequentially.a The two multiplexings are used separately with results shown with a box plot as in Fig.2.The fidelity achieved by random guessing is indicated by the dashed horizontal line.b The two multiplexings are used together.In all cases results of 100 random realizations of the reservoir, input and channel have been used.

Figure 4 .
Figure 4. Results for entangler task.For each different case 100 random networks were created, whereas the figure of merit is average logarithmic negativity.a Delay is varied for fixed N and all results are shown with a box plot as in Fig. 2. Delay 1 corresponds to nearest neighbors, 2 to next nearest neighbors and so on.b Median logarithmic negativity when N and delay are varied.See text for details.

Figure 5 .
Figure5.Partial generalizations.a Results for predictive quantum state preparation where the input is classical but target output is quantum.The reservoir is trained to prepare the quantum state according to future classical inputs which it must deduce from previous inputs.Specifically, the inputs are thermal states where the number of thermal excitations follows the well known Santa Fe time series, and the targets are squeezed vacuum states where the squeezing parameter is to coincide with the number of thermal excitations of a future input.b Results for von Neumann entropy detection where the input is quantum but the target output is classical.Here the reservoir is trained to estimate the determinants of input covariance matrices-which completely determines, e.g., the von Neumann entropy-for different delays.Random product states of M = 10 identical single mode states are used as input.If reservoir observables are available for a single timestep, then the von Neumann entropies of multiple previous inputs can be estimated with a very low error.