# Squeezing as a resource for time series processing in quantum reservoir
computing

Jorge García-Beni  Gian Luca Giorgi  Miguel C. Soriano  Roberta Zambrini
Instituto de Física Interdisciplinar y Sistemas Complejos (IFISC), UIB–CSIC
UIB Campus, Palma de Mallorca, E-07122, Spain

###### Abstract

Squeezing is known to be a quantum resource in many applications in metrology,
cryptography, and computing, being related to entanglement in multimode
settings. In this work, we address the effects of squeezing in neuromorphic
machine learning for time series processing. In particular, we consider a
loop-based photonic architecture for reservoir computing and address the
effect of squeezing in the reservoir, considering a Hamiltonian with both
active and passive coupling terms. Interestingly, squeezing can be either
detrimental or beneficial for quantum reservoir computing when moving from
ideal to realistic models, accounting for experimental noise. We demonstrate
that multimode squeezing enhances its accessible memory, which improves the
performance in several benchmark temporal tasks. The origin of this
improvement is traced back to the robustness of the reservoir to readout noise
as squeezing increases.

††preprint: APS/123-QED

##  I Introduction

Squeezing is a quantum phenomenon characterized by reduced light field
quadrature fluctuations below shot noise levels [1, 2, 3]. Initially employed
in fundamental quantum tests, such as Einstein–Podolsky–Rosen (EPR) paradox
experiments [4, 5], squeezing has emerged as a crucial resource in diverse
quantum technologies. Notably, squeezed states have been extensively utilized
in quantum metrology to enhance measurement sensitivity for parameter
estimation [6, 7], clock synchronization [8], and gravitational wave detection
[9]. Moreover, their role as a resource for quantum entanglement has been
harnessed for quantum cryptography protocols [10, 11]. In boson sampling
experiments, large multimode squeezed states have made it possible to achieve
a quantum advantage [12, 13]. Additionally, they serve as the primary resource
for universal measurement-based quantum computing in continuous variables (CV)
[14] through the generation of cluster states [15, 16, 17]. In the context of
quantum machine learning, squeezing is fundamental for CV quantum neural
networks to outperform their classical counterparts [18]. In this work, we
will focus on the favorable impact of squeezing on time-series prediction and
forecasting in the context of Quantum Reservoir Computing (QRC) [19].

Reservoir Computing (RC) constitutes an unconventional paradigm within the
realm of machine learning techniques rooted in recurrent neural networks [20,
21, 22]. Particularly tailored for time series processing, RC allows fast
learning with minimal training costs. The RC framework has demonstrated its
effectiveness in real-world scenarios including temporal prediction tasks [23,
24, 25, 26] as well as classification tasks [27, 28]. By harnessing the
information processing capabilities of high-dimensional dynamical systems, RC
concepts have seamlessly transitioned to physical substrates [29], with
photonic and optoelectronic implementations receiving attention for their
high-speed attributes [30, 31, 32, 33].

Recently, the scope of RC has expanded to encompass quantum systems,
capitalizing on their augmented Hilbert space for enhanced performance [19].
Notably, quantum enhancements in temporal tasks under ideal conditions have
been observed in both spin [34, 35] and photonic setups [36, 37]. Different
aspects influencing quantum reservoirs perfomance have been considered and the
effectiveness of complex task solving has been addressed considering different
evolution maps [34, 38, 39, 37], the role of statistics [40, 41], or different
quantum phases [42]. These investigations assumed ideal conditions,
attributing performance improvements to factors such as improved memory
properties, more favorable nonlinearities, and expanded accessible Hilbert
spaces.

However, QRC faces several challenges, with substantial attention directed
towards addressing the presence of noise in output observables [43, 44, 45,
46, 47]. While readout noise is relevant to classical RC as well [48, 49], it
acquires heightened significance in QRC due to intrinsic sampling noise
arising from the stochastic nature of quantum measurements. This noise
significantly hinders potential quantum enhancements [44, 45, 46]. The
strategy of monitoring the output accounting on the effects of quantum
measurement for temporal tasks is also particularly critical [50, 45].
Preserving quantum advantage within these non-ideal circumstances is needed to
ensure the viability of QRC protocols.

In Ref. [50] it was shown how weak –instead of projective– measurements allow
continuous monitoring in QRC –instead of buffering inputs and rewind the
reservoir dynamics. Continuous homodyne monitoring in QRC was addressed in
Ref. [45] in a photonic platform. Here, a loop-based photonic platform was
proposed, suitable for online time series processing. These works show the
possibility of sustaining the RC performance in the quantum setting in the
presence of non-ideal measurement conditions. Motivated by these results, in
this work, we adopt a slightly simplified version of the photonic platform
used in Ref. [45] and investigate how the presence of squeezing in the optical
cavity can improve the performance of the reservoir. We will numerically
analyze the role of squeezing, addressing active and passive coupling terms in
a photonic network, in both linear and nonlinear tasks. An analytical argument
is made to justify such an improvement.

The paper is structured as follows: in Sect. II, the general framework of RC
as well as a detailed description of our photonic platform are exposed. In
Sect. III the simulation results for some benchmark RC tasks under non-ideal
noise conditions are shown: we test the linear memory of the system (Sect.
III.1) as well as the non-linear memory (Sect. III.2) and its performance on
time-series forecasting (III.3). In Sect. IV we show numerical evidence to
explain the noise robustness improvement caused by squeezing in the previous
tasks. Finally, conclusions are given in Sect. V.

##  II Loop-based architecture

###  II.1 Reservoir computing

RC architectures are mainly composed of three distinct layers: the input
layer, the reservoir, and the readout layer [51]. The external signal is
encoded and fed into the system in the input layer. This is done sequentially
at each time step. The reservoir layer (or just reservoir) is usually a
complex dynamical system that applies a non-linear map to the inputs. The
reservoir must retain short-term memory of previous inputs to be able to
perform temporal tasks. This short-term or fading memory, together with the
echo state property, is part of the universality proofs of RC [52]. The
readout layer is then made of a certain number of reservoir observables, which
are monitored sequentially after the input is encoded. The output from this
layer is a linear combination of the measured observables. Supervised learning
is performed by optimizing this linear combination of the output to yield the
desired target.

In more detail, if we have a training set consisting of a sequence of $L$
inputs $\left\\{s_{1},s_{2},\dots,s_{L}\right\\}$, at time step $k$ the input
$s_{k}$ is encoded and introduced into the reservoir. If we call
$\mathbf{x}_{k}$ to the reservoir degrees of freedom, we can write the
reservoir map at time step $k$ as

$\mathbf{x}_{k}=\mathcal{H}\left(\mathbf{x}_{k-1},s_{k}\right).$ (1)

This map is fixed throughout the whole protocol. For the readout layer at time
step $k$, we use $\mathbf{O}_{k}$ as readout observables (functions of
$\mathbf{x}_{k}$). The reproduced function at each time is obtained by
performing a linear regression on the readout observables,

$y_{k}=w_{0}+\mathbf{O}_{k}^{\top}\mathbf{w}\\\
=\left(1,\mathbf{O}_{k}^{\top}\right)\left(\begin{array}[]{c}w_{0}\\\
\mathbf{w}\end{array}\right),$ (2)

where the weight vector
$\mathbf{W}=\left(w_{0},\mathbf{w}^{\top}\right)^{\top}$ is optimized through
training examples. The way training works is as follows: for the sequence of
$L$ training inputs, we define the matrices

$V=\left(\begin{array}[]{cc}1&\mathbf{O}_{1}^{\top}\\\
1&\mathbf{O}_{2}^{\top}\\\ \vdots&\vdots\\\
1&\mathbf{O}_{L}^{\top}\end{array}\right);\quad\mathbf{y}=\left(\begin{array}[]{c}y_{1}\\\
y_{2}\\\ \vdots\\\ y_{L}\end{array}\right),$ (3)

so Eq. y-linear-reg can be rewritten as $\mathbf{y}=V\mathbf{W}$. If we want
the output $\mathbf{y}$ to get as close as possible to the desired target
function, $\bar{\mathbf{y}}$, we choose the set of weights in order to
minimize the mean square error (MSE),

$\text{MSE}\left(\mathbf{y},\bar{\mathbf{y}}\right)=\frac{1}{L}\sum_{k=1}^{L}\left(y_{k}-\bar{y}_{k}\right)^{2}.$
(4)

It can be shown that the optimal set of weights to reach the minimum of Eq.
mse-eq are the ones obtained through
$\mathbf{W}_{\text{opt}}=V^{\text{MP}}\bar{\mathbf{y}}$, where
$V^{MP}=\left(V^{\top}V\right)^{-1}V^{\top}$ is the Moore-Penrose inverse of
$V$ [22]. The higher the value of $L$ the more precise our estimation of the
optimal weights will usually be. Once the system has been trained, we consider
a (smaller) test set of $L^{\prime}$ inputs to be fed into the reservoir
afterwards. The RC performance is then checked with this test set of new
unseen data using the MSE metric from Eq. mse-eq.

###  II.2 Description of the platform

Our architecture works in the CV quantum optical regime. The physical
substrate is an $N$-mode optical pulse traveling through a closed optical loop
or cavity ($N$ denotes the size of the reservoir). The $N$-mode internal
degrees of freedom inside each pulse can be attained via, for example,
frequency multiplexing [53, 54, 16, 55, 56]. In optics, frequency multiplexing
has already been shown to be a useful strategy for classical RC [57]. In our
approach, the external information is injected from a pulse-generating light
source, which provides squeezed vacuum states. The input sequence is encoded
in the squeezing phases of the input pulses (one input value for each pulse),
depicted as source in Fig. 1. Fast, accurate, and reconfigurable phase-setting
devices have already been used in experiments with great impact [13]. Each
input pulse is coupled to the loop pulse using a beam-splitter (BS), with
reflectivity $R$ shown in Fig. 1, yielding two output pulses. One of them
remains in the loop and gives feedback to the next iteration (creating a
quantum memory). In this way, the reservoir can retain information from
previous inputs without the need for external memory. The remaining output
pulse is passed to a detector that measures each mode and uses the obtained
observables for the readout layer. The fraction of light that remains in the
cavity on each round trip is determined by the BS reflectivity $R$.

Inside the cavity, a nonlinear medium is placed (NL in Fig. 1), which applies
a dynamical transformation to the loop pulse each time it passes through. This
creates a complex optical network [54, 58, 59, 55, 60] and can be modeled by a
Hamiltonian that is quadratic in the field operators,

$\hat{H}=\frac{1}{2}\sum_{i,j=1}^{N}\left(\alpha_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}+\beta_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}+\text{h.c.}\right),$
(5)

where $\hat{a}_{i}$ ($\hat{a}_{i}^{\dagger}$) is the annihilation (creation)
operator of mode $i$. The coupling terms $\alpha_{ij}$ and $\beta_{ij}$ encode
different network topologies and lead to entanglement among modes inside the
loop pulse, at each round trip. If all the terms $\beta_{ij}=0$, the dynamical
transformation is called passive, whereas if there are any $\beta_{ij}\neq 0$
the transformation is active. Active transformations in CV quantum optics do
not conserve the average number of photons of the quantum state and are known
to generate squeezing, the main resource for entanglement [2, 3]. It is
important to note that a passive cavity also produces entangled states because
the external input pulses are already squeezed (even though it does not
generate additional squeezing).

Figure 1: Scheme of the loop-based architecture.

The detector from Fig. 1 performs homodyne measurements to every mode in the
incoming pulses. Concretely, the measured operator are the $x$-quadrature
$\hat{x}_{i}=\frac{1}{\sqrt{2}}\left(\hat{a}_{i}+\hat{a}_{i}^{\dagger}\right)$
($i=1,2,\dots,N$). From these measurements, the moments of the operator vector
$\hat{\mathbf{x}}=\left(\hat{x}_{1},\dots,\hat{x}_{N}\right)^{\top}$ can be
computed and used as observables for the readout layer. As we are injecting
squeezed vacuum states, only even-order moments are considered. For the tasks
shown in Sect. III, the chosen set of observables is composed of second and
fourth-order moments. Concretely, the chosen set is
$\left\\{\left\langle\hat{x}_{i}\hat{x}_{j}\right\rangle,\left\langle\hat{x}_{i}^{2}\hat{x}_{j}^{2}\right\rangle,\left\langle\hat{x}_{i}^{3}\hat{x}_{j}\right\rangle\right\\}_{i,j=1}^{N}\equiv\left\\{O_{l}\right\\}_{l=1}^{N(3N+1)/2}$,
which has a total of $N(3N+1)/2$ observables. In Gaussian states, fourth-order
moments can be written as nonlinear functions of the second-order moments. In
our case, they are useful for enhancing the accessible nonlinear terms, which
are relevant for several tasks and come at no experimental expenditure [45].
We note that the readout size scales quadratically with the number of modes
$N$. Usually, in optical reservoir computing the dimensionality scales
linearly with $N$, as the information is encoded in field amplitudes [32, 33,
61]. By introducing the inputs in the field quantum fluctuations we access a
broader dimensional space by generating mode correlations and entanglement,
which allows more complex information processing in relatively smaller
reservoirs [36, 45]. To access the averaged values of the field correlations
we consider averaging over an ensemble of realizations [34].

Each RC time step encompasses the whole process we have detailed: BS coupling
of the input and loop pulses, detection of the output pulse, and
transformation of the loop pulse under the non-linear crystal. There are as
many time steps as samples in the input sequence, and they will be labeled
with the letter $k$. So at the $k$-th time step, the input $s_{k}$ will be
encoded and introduced into the system and the vector of observables
$\mathbf{O}^{(k)}=\left(O_{1}^{(k)},O_{2}^{(k)},\dots\right)$ is measured and
used for the readout layer.

The platform is based on the proposal in Ref. [45], where real-time
information processing was reported. The main novelty here is that we
specifically address the role of a squeezing reservoir, with both active and
passive transformation. The goal is to assess the importance of quantum
resources for the performance of QRC. In order to simplify the experimental
footprint consider a single NL crystal. Indeed this has not significant effect
when processing past inputs. This design can be also adapted for single-loop
ensemble processing adding a fiber, as shown in [45].

##  III Reservoir computing tasks under additive noise

Noise in the readout layer is known to be significantly detrimental to RC
performance [48, 49]. Some strategies have been developed in different
architectures to make reservoirs more robust to noise [48, 43]. In this
section, we study the effect of the amount of squeezing produced by the cavity
crystal on the noise robustness of the platform in the readout layer and
compare it to the one obtained by tuning the BS reflectivity. In our
simulations, readout noise is included as additive fluctuations in the
measured observables, $\mathbf{O}_{\text{meas}}^{(k)}$, concretely

$\mathbf{O}_{\text{meas}}^{(k)}=\mathbf{O}_{\text{ideal}}^{(k)}+\mathbf{\mathcal{E}}^{(k)}(0,\sigma_{\text{noise}}^{2}),$
(6)

where $\mathbf{O}_{\text{ideal}}^{(k)}$ stands for the observable that would
be measured in the ideal case of zero fluctuations, which would correspond to
an infinite number of measurements, and $\mathbf{\mathcal{E}}^{(k)}$ stands
for the additive noise vector. For the noise, we model it as normally
distributed fluctuations of variance equal to $\sigma_{\text{noise}}^{2}$
applied to the measured quadratures. Although the added noise is absolute (it
does not depend on the magnitude of the observables), we can use the vacuum
noise variance ($\sigma_{\text{vacuum}}^{2}=1/2$ in our case) or shot noise as
a relative measure of the additive noise intensity. In that regard, added
noise of variance equal to $0.1$ would be equivalent to $20\%$ of vacuum
fluctuations.

In our simulations, we generate every crystal Hamiltonian, Eq. general-hamil,
randomly with the condition that every one of its supermodes is squeezed by
$e^{-r}$ (see App. A for details). In every realization, the modes of the
input pulses are squeezed with a fixed squeezing strength,
$r_{\text{input}}=2$ (approximately 8.7 dB). The encoding function of the
squeezing phase, $\phi_{k}=f(s_{k})$, is tuned depending on the task we are
considering. Concretely, we consider the family of linear functions
$\phi_{k}^{(m)}=m\pi s_{k}$. In this respect, the smaller the value of $m$,
the better the reservoir is at reproducing linear and quadratic functions of
$s_{k}$. For increasing values of $m$, higher non-linear contributions become
more relevant [45, 36].

We consider three temporal tasks: the linear memory task, the nonlinear
autoregressive moving average (NARMA) task, and the forecasting of the Mackey-
Glass chaotic time series [62, 63]. These three tasks provide a broad overall
picture of the properties of the reservoir memory for both linear and
nonlinear computations. After applying the training protocol described in Sect
II.1 for a given target function, we check the performance of the trained
reservoirs on an additional test input sequence. For evaluation, we mainly use
the mean-square error, Eq. mse-eq, after training optimization,
$\text{min}_{\mathbf{W}}\text{MSE}\left(\mathbf{y},\bar{\mathbf{y}}\right)$.
To get the MSE to be normalized between 0 and 1, we normalize both the target
and the reproduced function data to zero mean and one standard deviation.

###  III.1 Linear memory task

The linear memory task is the simplest way to check the accessible memory of
our reservoir in the presence of noise. We train the reservoir to reproduce
past entries from the input series. That is, we consider the target function

$\bar{y}_{k}(\mathbf{s},d)=s_{k-d},$ (7)

where we aim to to reproduce, at time step $k$, the input that was introduced
$d$ time steps in the past (or the input at delay $d$). For simulations, we
consider an input sequence composed of random entries from a uniform
distribution in the interval $[-1,1]$. We tune the squeezing angle encoding to
be $\phi_{k}=\pi s_{k}/4$ to maximize linear contributions. For visualization
purposes, we use the linear capacity, defined as
$C(d)=1-\text{min}_{\mathbf{W}}\text{MSE}\left(\mathbf{y},\bar{\mathbf{y}}(d)\right)$,
to test the performance of this task. The vector $\bar{\mathbf{y}}(d)$ is
composed of the target function from Eq. linear-target-function at each time
step as a function of the delay $d$. Training and test set sizes are set to be
4000 and 1000, respectively.

Figure 2: Linear memory under additive noise: (a-b) linear capacity as a
function of the delay for different values of squeezing (green for $r=0$,
orange for $r=0.75$ and purple for $r=1.5$) and different reflectivities:
$R=0.75$ in (a) and $R=0.9$ in (b). In both figures, the noise variance is
$\sigma_{\text{noise}}^{2}=10^{-2}$. The curves are taken from averaging among
100 different random realizations and the shadows depict the standard
deviation. (c) delay at which the linear capacity drops below $0.9$ as a
function of $\sigma_{\text{noise}}^{2}$, including the ideal case
($\sigma_{\text{noise}}^{2}=0$) for different values of squeezing (different
colors) and different values of the reflectivity (line format).

Figures 2a and 2b show the linear capacity as a function of the delay $d$ on
the target function from Eq. linear-target-function for different values of
the reflectivity ($R=0.75$ in Fig. 2a and $R=0.9$ in Fig. 2b). In both plots,
the noise variance is $10^{-2}$ ($2\%$ of vacuum fluctuations). We see that
for both reflectivities, increasing the cavity squeezing improves the
attainable memory. Having a higher reflectivity provides a longer ‘tail’ in
the linear capacity at the expense of reducing it for small and intermediate
delays. This is the effect of the smaller amount of light leaving the cavity
for increasing values of $R$. In Fig. 2c the delay at which the linear
capacity drops below $0.9$ (which we also call the delay cut) is plotted as a
function of the noise variance. We find that cavity squeezing provides
significant noise robustness. For the reflectivity $R=0.75$, adding cavity
squeezing equal to $r=1.5$ provides a high linear capacity beyond delay $10$
for a noise intensity of $\sigma_{\text{noise}}^{2}=0.1$. In Fig. 2c the
drawback of increasing the BS reflectivity can also be noted: for small noise
intensity ($\sigma_{\text{noise}}^{2}=10^{-3}$) and no cavity squeezing (light
color), having $R=0.9$ (dashed line) improves the delay cut compared to
$R=0.75$ (solid line), as increasing the reflectivity improves the memory
robustness to noise. However, when the noise intensity is increased, the delay
cut for $R=0.9$ drops below the one for $R=0.75$. This is because the higher
the reflectivity, the less light leaves the cavity and travels to the
detector. If the readout fluctuations are large compared to the intensity of
the light coming from the loop, the accessible memory of the reservoir will be
severely degraded. The cavity squeezing increases significantly the accessible
linear memory in the presence of a large noise intensity.

###  III.2 NARMA10 task

In this section, we are analyzing the performance of the reservoir for the
NARMA task, which requires high linear and non-linear memory. It is one of the
most common benchmark tasks and has been used to test several QRC proposals
[42, 43]. In this article, we consider the NARMA10 task, with a target
function at time $k$ being

$\displaystyle\bar{y}_{k}$ $\displaystyle=$
$\displaystyle\alpha\bar{y}_{k-1}+\beta\bar{y}_{k-1}\sum_{i=1}^{10}\bar{y}_{k-i}+\gamma
u_{k-1}u_{k-10}+\delta,$ (8) $\displaystyle u_{k}$ $\displaystyle=$
$\displaystyle\mu+\nu s_{k},$ (9)

where the default constant parameters are set to
$(\alpha,\beta,\gamma,\delta)=(0.3,0.05,1.5,0.1)$. The function input
parameters $(\mu,\nu)$ are chosen to be $(0,0.2)$, where $s_{k}$ are taken
from a uniform random distribution in the interval $[-1,1]$ (they are the
source inputs). For the reservoir to be able to have a good performance on the
NARMA10 task it is necessary for it to be able to have a high linear capacity
up to delay 10 and a low error reproducing the function $y_{k}=s_{k-1}\cdot
s_{k-10}$ [64].

To perform this task we consider the same input encoding as in Sect. III.1,
namely $\phi_{k}=\pi s_{k}/4$. To test the performance we use the MSE defined
in Eq. mse-eq. In Figs. 3a and 3b the performance is shown, comparing the
effect of cavity squeezing (x-axis in Fig. 3a) and BS reflectivity (x-axis in
Fig. 3b). Three different noise scenarios are considered in both figures: the
ideal case (blue boxes) and noise variance equal to $10^{-2}$ (green boxes)
and $10^{-1}$ (pink boxes). In the ideal case, the optimal values of cavity
squeezing and BS reflectivity are found to be $r=0$ and $R=0.5$. It can be
seen that increasing either cavity squeezing or reflectivity degrades the
performance in the absence of noise. The reason for this is that, for the
chosen encoding and observables, increasing these two parameters increases
linear memory at the expense of quadratic memory, which is also very relevant
for the NARMA task. This balance among linear and non-linear memory is well
known in the field of RC [65, 66, 36].

Moving to the realistic case of a finite number of measurements, for a noise
intensity of $\sigma_{\text{noise}}^{2}=10^{-2}$, the improvement of the
cavity squeezing over the reflectivity is clearly seen. In that noise
scenario, increasing both parameters improves the performance up to an optimal
value ($r\sim 1$ and $R\sim 0.8$). However, in the passive cavity case, the
optimal value has a higher error in comparison to the active cavity and is
further away from the ideal case error. For higher values of the noise,
$\sigma_{\text{noise}}^{2}=10^{-1}$, most examples completely fail at
attempting to reproduce the NARMA function. Only when the cavity squeezing
$r\gtrsim 1$, the MSE drops below 1. This is due to the fact that in most
scenarios the noise level makes it impossible for the reservoir to resolve
inputs with a delay of 10 or more as its intensity becomes comparable with the
value of the system observables, giving a bad signal-to-noise ratio [45].

Even though there is a counterbalance between linear and quadratic memory, and
thus increasing the cavity squeezing and the BS reflectivity is detrimental to
the performance in ideal scenarios, in the presence of readout noise the
active cavities with high squeezing outperform the rest. Indeed, the only case
where the performance of the noisy reservoir comes close to the ideal case is
when the cavity squeezing is higher than $r=1$. This provides an interesting
example where the role of a given resource needs to be addressed beyond ideal
settings, as benefits could arise in the presence of noise.

Figure 3: Performance on NARMA10 task: box plot of the mean square error (MSE)
of the NARMA10 task as a function of (a) squeezing and (b) reflectivity for
different values of the noise variance (different colors). For Fig. 3a we take
the reflectivity equal to $R=0.5$ and in Fig. 3b the cavity squeezing is equal
to zero. For a given value of the x-axis, the boxes for each noise scenario
are split to avoid overlapping.

###  III.3 Time series prediction of a chaotic signal

One of the main applications of RC is time series forecasting, and thus in
this section, we consider the task of forecasting the Mackey-Glass time series
[67]. The differential equation that describes the dynamics of the signal is

$\dot{s}(t)=-0.1s(t)+\frac{0.2s(t-\tau)}{1+s(t-\tau)^{10}},$ (10)

where for $\tau=17$ the time series is chaotic [62, 63]. For the input
sequence, we have sampled the solutions to Eq. mackey-glass-diffeq with time
resolution $t_{r}=3$ [61], so that $s_{k}=s(t_{0}+kt_{r})$ (the initial
conditions are chosen randomly). The target function for the training will be
to predict the next input in the sequence, that is, $\bar{y}_{k}=s_{k+1}$. For
this specific task, we use the input encoding $\phi_{k}=\pi s_{k}$, which
provides higher nonlinear memory [45]. Once the reservoir has been trained to
predict the next value of the signal, we can feed the predicted values as new
inputs for the protocol. The reservoir thus, ideally, faithfully reproduces
the chaotic signal without the need for new input data. We call this protocol
autonomous driving.

In Fig. 4a the autonomously driven signal evolution is plotted for two
different values of the cavity squeezing (green curve for $r=0$ and blue curve
for $r=1.25$) while the noise variance is kept at
$\sigma_{\text{noise}}^{2}=10^{-1}$. The real signal is shown as a black curve
for comparison. We see that the active cavity performs a much better
prediction than the passive one in the long term, achieving accurate
predictions up to 50 time steps. The passive cavity reservoir cannot overcome
the effects of noise, and thus the prediction performance drops dramatically.
In Figs. 4b and 4c, we compare the Mackey-Glass chaotic attractor, black dots,
to $3000$ values from single realizations of trained reservoirs with a passive
cavity (Fig. 4b, green dots) and an active cavity with $r=1.25$ (blue dots).
Also in this case, we see that the reservoir with cavity squeezing is able to
approximately reproduce the attractor, while the passive reservoir is not.

Figure 4: Mackey-Glass time series prediction under additive noise: (a) time
series as a function of the reservoir time steps: the black curve shows the
real-time series while the green and blue curves show the autonomously driven
predictions for $r=0$ and $r=1.25$, respectively (averages taken among 100
realizations, shadows depicting the standard deviation). (b-c) chaotic
attractor in the phase space $y(t)$-$y(t-6)$: the black points provide the
real attractor while the green and blue dots provide the results of an
autonomously driven realization with $r=0$ and $r=1.25$, respectively. In
every simulation, the BS reflectivity is set to $R=0.75$ and the noise
variance is $\sigma_{\text{noise}}^{2}=10^{-1}$.

##  IV Accessible memory enhancement

In this section, we will explain in detail the reason behind the performance
improvement due to the cavity squeezing that we have shown in the previous
sections. From a physical point of view, the BS and the cavity crystal can
induce competing effects on the reservoir memory. The BS causes a loss of
photons in the loop pulse. On the other hand, the non-linear processes that
arise from the interactions inside the crystal can increase the energy of the
loop pulse. Concretely, active crystals (those that produce squeezing)
increase the total photon number inside the pulse, as opposed to passive
crystals, which maintain it constant. This energy enhancement from active
crystals counteracts BS losses and helps retain information inside the loop
pulse for longer times.

To quantify how these effects contribute constructively to the functioning of
the QRC, we study the time evolution of the loop pulse during the protocol. At
each round trip, the field quadratures of the loop pulse transform via the
symplectic matrix $A=\sqrt{R}S$, where $R$ stands for the BS reflectivity and
$S$ is the symplectic matrix modeling the evolution of the pulse inside the
cavity crystal (see Apps. A and B for details on the matrices $S$ and $A$,
respectively). The memory retention of our reservoir is directly related to
the powers of the symplectic matrix $A$ (App. B.1) and can be quantified using
the spectral norm of $A^{d}$ (written as $\lVert A^{d}\rVert_{2}$), in which
$d$ resembles the delay of the input information (Eq. spectral-norm from App.
B.2). In that regard, the faster $\lVert A^{d}\rVert_{2}$ decays to zero, the
smaller the memory retention of our reservoir will be (and vice-versa). It can
be analytically shown that, if the crystal inside the cavity is passive, then
$\lVert A^{d}\rVert_{2}=R^{d/2}$ and, if it is active, then $\lVert
A^{d}\rVert_{2}\geq R^{d/2}$. From these equations, we infer that active
cavity crystals can improve the memory retention of the reservoir (similarly
to the direct effect of the BS reflectivity). In Fig. 5 we show the values of
$\lVert A^{d}\rVert_{2}$ as a function of the delay for the randomly generated
Hamiltonians that we considered in Sect. III. In Fig. 5a the spectral norm of
$A$ is plotted for active cavities with different values of cavity squeezing,
$r$, while in Fig. 5b the same is shown for passive cavities with different
values of the BS reflectivity. We see how the negative slope of the spectral
norm decreases as we increase either the cavity squeezing (Fig. 5a) or the BS
reflectivity (Fig. 5b). This means that the magnitude of the delayed input
information decays more slowly, and thus it may be reproduced by our trained
reservoir more easily. In the presence of readout noise, input information
decaying more slowly with time translates into more accessible memory and
improves the reservoir’s robustness to noise This is the reason why cavity
squeezing improves the performance in every benchmark task shown in Sect. III.
While both $R$ and $A$ contribute to such a memory enhancement, increasing the
BS reflectivity has the drawback of reducing the cavity light reaching the
detector (lowering the signal-to-noise ratio) and thus is not as effective for
increasing the noise robustness.

Figure 5: Spectral norm of loop matrix dynamics: spectral norm of the loop
dynamical matrix $A$ to the power of the delay $d$ for (a) active
transformations with $R=0.5$ and different values of the squeezing, $r$, and
(b) passive transformations ($r=0$) for different values of the reflectivity.
Each curve shows the median from 100 different realizations while the shades
range from the first to the ninth decile.

##  V Conclusion

In light of the significant achievements in classical RC [68], photonic
platforms are emerging as promising candidates for quantum implementations.
Commendable features are, for instance, fast processing rates and low
decoherence at room temperature [19]. Different photonic QRC platforms have
already been theoretically explored and show improvements due to the enlarged
Hilbert space [36, 69] as well as the ability of real-time processing without
the use of external memories [45].

Detrimental effects of readout noise have been discussed both in classical RC
[48, 49] and in QRC settings [43, 44, 45]. In the case of quantum reservoirs,
the problem is even more profound, as readout noise is theoretically
unavoidable due to the stochastic nature of quantum measurements, which
produce statistical fluctuations that can hinder any possible quantum
advantage [19, 44, 45]. In this paper, we demonstrate the performance-
enhancing potential of quantum squeezing (active cavity) applied to a vacuum
state quantum memory (loop pulse) to overcome noise in realistic scenarios.
This establishes squeezing as a quantum resource for accessing the enlarged
space of quantum correlations and entanglement and for improving the
performance in relevant benchmark tasks, either predictive or requiring
memory. Even though tuning the BS reflectivity also improves memory retention,
as in [45], increasing cavity squeezing is shown to be a preferred method to
improve the reservoir robustness under adverse noise conditions.
Interestingly, the effect of squeezing on the QRC performance when accounting
for measurement noise completely deviates from predictions under ideal
conditions.

In summary, state-of-the-art frequency multiplexed quantum networks [16, 55,
56, 60] represent a powerful set-up for near-term experimental implementations
of QRC. Our results serve as a guide for the experimental design laying the
foundations for photonic QRC in CV in realistic noisy scenarios and exploiting
quantum resources.

###### Acknowledgements.

We acknowledge the Spanish State Research Agency through the QUARESC project
(PID2019-109094GB-C21 and -C22/ AEI / 10.13039/501100011033) and through the
Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence
in R&D (CEX2021-001164-M funded by the MCIN/AEI/10.13039/501100011033) and. We
also acknowledge funding by CAIB through the QUAREC project (PRD2018/47).
J.G-B. is funded by the Conselleria d’Educació, Universitat i Recerca of the
Government of the Balearic Islands with grant code FPI/036/2020. G.L.G. is
funded by the Spanish MEFP/MiU and co-funded by the University of the Balearic
Islands through the Beatriz Galindo program (BG20/00085). This work has been
financially supported by MINECO through the QUANTUM ENIA project, and by EU
through the RTRP - NextGenerationEU within the framework of the Digital Spain
2026 Agenda. The CSIC Interdisciplinary Thematic Platform (PTI) on Quantum
Technologies in Spain is also acknowledged.

##  Appendix A Symplectic formalism of the crystal dynamics

A CV quantum state is completely determined, in the Heisenberg picture, by the
statistics of the quadrature vector
$\hat{\mathbf{Q}}=\left(\hat{x}_{1},\hat{p}_{1},\dots,\hat{x}_{N},\hat{p}_{N}\right)^{\top}$,
where
$\hat{x}_{i}=\frac{1}{\sqrt{2}}\left(\hat{a}_{i}^{\dagger}+\hat{a}_{i}\right)$
and
$\hat{p}_{i}=\frac{i}{\sqrt{2}}\left(\hat{a}_{i}^{\dagger}-\hat{a}_{i}\right)$
($i=1,\dots,N$) are, respectively, the amplitude and phase quadratures of each
mode. For quadratic Hamiltonians as in Eq. general-hamil the evolution of the
quadrature vector can be written as

$e^{i\hat{H}t}\hat{\mathbf{Q}}e^{-i\hat{H}t}=S_{\hat{H}}(t)\hat{\mathbf{Q}},$
(11)

where $S_{\hat{H}}$ is a $2N\times 2N$ symplectic matrix [2, 3]. For the sake
of clarity, we will drop the $\hat{H}$ subscript (and the time dependency) and
call the symplectic matrix simply as $S$. Every symplectic matrix admits a
Bloch-Messiah decomposition [70, 71],

$\displaystyle S$ $\displaystyle=$ $\displaystyle U\Delta_{\mathbf{\xi}}V,$
(12) $\displaystyle\Delta_{\mathbf{\xi}}$ $\displaystyle=$
$\displaystyle\text{diag}\left(\xi_{1},\xi_{1}^{-1},\xi_{2},\xi_{2}^{-1},\dots,\xi_{N},\xi_{N}^{-1}\right),$
(13)

where $U$ and $V$ are orthogonal symplectic matrices and $\xi_{i}^{-1}$
($\xi_{i}$) provides the squeezing (anti-squeezing) applied to the $i$-th
supermode (we consider $\xi_{i}\geq 1\ \forall i$). Thus the squeezing and
anti-squeezing values obtained from the Bloch-Messiah decomposition correspond
to the singular values of $S$. So the transformation is passive iff all the
singular values of $S$ are equal to 1 ($\Delta_{\mathbf{\xi}}$ is equal to the
identity). It is trivial to show that if $S$ is passive, then all its powers
$S^{d}$ are also passive.

For the simulations we have performed in this manuscript, we have considered
the non-linear crystal squeezes every supermode by the same amount, so
$\xi_{i}^{-1}\equiv e^{-r/2}$ ($\forall i$), where $r$ denotes the squeezing
strength per mode applied by the Hamiltonian. So for our simulations we set
$\Delta_{\xi}=\bigoplus_{i=1}^{N}\text{diag}\left(e^{r/2},e^{-r/2}\right)$ and
choose the matrices $U$ and $V$ randomly to generate the symplectic matrix
$S$, applying Eq. bloch-messiah-S.

##  Appendix B QRC dynamics

###  B.1 Recursive equations

In this section, we will work out the equations describing the reservoir
dynamics to gain insight into the effect of cavity squeezing. We start by
considering the quadrature operators of the pulses after the BS coupling at
time step $k$. We write them as

$\displaystyle\hat{\mathbf{Q}}_{\text{out}}^{(k)}$ $\displaystyle=$
$\displaystyle\sqrt{1-R}\hat{\mathbf{Q}}_{\text{loop}}^{(k)}-\sqrt{R}\hat{\mathbf{Q}}_{\text{in}}^{(k)},$
(14) $\displaystyle\hat{\mathbf{Q}}_{\text{loop}}^{(k)\prime}$
$\displaystyle=$ $\displaystyle
S\left(\sqrt{R}\hat{\mathbf{Q}}_{\text{loop}}^{(k)}+\sqrt{1-R}\hat{\mathbf{Q}}_{\text{in}}^{(k)}\right),$
(15)

where $\hat{\mathbf{Q}}_{\text{in}}$ and $\hat{\mathbf{Q}}_{\text{loop}}$ are
the quadrature operator vector of the input and the cavity pulse
(respectively) before the coupling, $R$ is the BS reflectivity and $S$ is the
symplectic transformation that the nonlinear medium applies. To find an
expression for the observables we compute the expressions of the covariance
matrix,
$\Gamma_{\text{loop}}^{(k)}=\left\langle\hat{\mathbf{Q}}_{\text{loop}}^{(k)}\hat{\mathbf{Q}}_{\text{loop}}^{(k)\top}\right\rangle-\left\langle\hat{\mathbf{Q}}_{\text{loop}}^{(k)}\right\rangle\left\langle\hat{\mathbf{Q}}_{\text{loop}}^{(k)}\right\rangle^{\top}$,
where the averages are taken from an ensemble of realizations. As we are
working in squeezed vacuum, the covariance matrix and the second-order moments
are identical. Also, as we are dealing with Gaussian states, the higher-order
moments are functions of the first and second-order moments of the
quadratures. It can be shown that the covariance matrix of the cavity pulse at
the next time step can be written as

$\Gamma_{\text{loop}}^{(k+1)}=RS\Gamma_{\text{loop}}^{(k)}S^{\top}+(1-R)S\Gamma_{\text{in}}^{(k)}S^{\top},$
(16)

which can be further expanded by recursion to yield

$\Gamma_{\text{loop}}^{(k+1)}=(1-R)S\left[\sum_{d=0}^{\infty}A^{d}\Gamma_{\text{in}}^{(k-d)}\left(A^{d}\right)^{\top}\right]S^{\top},$
(17)

for $A=\sqrt{R}S$. The details of the derivations of Eqs. CM-loop-simple and
CM-loop-expanded can be found in our previous work [45]. In Eq. CM-loop-
expanded, the index $d$ resembles the delay of the input.

###  B.2 Input memory decay

From Eq. CM-loop-expanded we can see that the magnitude decay of the input
information ($\Gamma_{\text{in}}$) depends on the powers of $A$, $A^{d}$,
where $d$ denotes the delay of the encoded input. The accessible memory of the
reservoir is determined by how much this input information decays each round
trip. To quantify this decay we use the spectral norm or L2-norm, defined as

$\lVert A\rVert_{2}=\max_{|\mathbf{x}|\neq
0}\frac{|A\mathbf{x}|}{|\mathbf{x}|}\equiv\xi_{\text{max}}(A),$ (18)

where $|\cdot|$ stands for the usual euclidian norm of a vector and
$\xi_{\text{max}}(A)$ is the maximum singular eigenvalue of $A$ [72]. From the
definition of $A$, we have that $\lVert A^{d}\rVert_{2}=R^{d/2}\lVert
S^{d}\rVert_{2}$. If the cavity is passive, we have that $\lVert
S^{d}\rVert_{2}=\lVert S\rVert_{2}=1$, and so $\lVert
A^{d}\rVert_{2}=R^{d/2}$. If the cavity is active we have that $\lVert
A^{d}\rVert_{2}\geq R^{d/2}$. For the RC protocol to work, it is important
that both the echo state property and the fading memory condition hold [73].
It can be shown that these conditions are fulfilled iff $\rho(A)<1$, where
$\rho(\cdot)$ stands for the spectral radius of a matrix [36]. Moreover, if
$\rho(A)<1$ then the following is also true:

$\lim_{d\to\infty}\lVert A^{d}\rVert_{2}=0.$ (19)

This ensures that the energy of the system does not diverge. We checked that
the condition $\rho(A)<1$ is held in every simulation in the manuscript.

## References

  * [1] Alessandro Ferraro, Stefano Olivares, Matteo GA Paris, et al.  Gaussian states in quantum information.  Bibliopolis, Napoli, 2005.
  * [2] Gerardo Adesso, Sammy Ragy, and Antony R. Lee.  Continuous variable quantum information: Gaussian states and beyond.  Open Systems & Information Dynamics, 21(01n02):1440001, Mar 2014\.
  * [3] Alessio Serafini.  Quantum continuous variables: a primer of theoretical methods.  CRC Press, 2017.
  * [4] M. D. Reid and P. D. Drummond.  Quantum correlations of phase in nondegenerate parametric oscillation.  Phys. Rev. Lett., 60:2731–2733, Jun 1988.
  * [5] Z. Y. Ou, S. F. Pereira, H. J. Kimble, and K. C. Peng.  Realization of the einstein-podolsky-rosen paradox for continuous variables.  Phys. Rev. Lett., 68:3663–3666, Jun 1992.
  * [6] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Quantum-enhanced measurements: Beating the standard quantum limit.  Science, 306(5700):1330–1336, Nov 2004.
  * [7] Adriano A. Berni, Tobias Gehring, Bo M. Nielsen, Vitus Händchen, Matteo G. A. Paris, and Ulrik L. Andersen.  Ab initio quantum-enhanced optical phase estimation using real-time feedback control.  Nature Photonics, 9(9):577–581, Sep 2015.
  * [8] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Quantum-enhanced positioning and clock synchronization.  Nature, 412(6845):417–419, Jul 2001.
  * [9] J. Aasi, J. Abadie, B. P. Abbott, R. Abbott, et al.  Enhanced sensitivity of the ligo gravitational wave detector by using squeezed states of light.  Nature Photonics, 7(8):613–619, Aug 2013.
  * [10] Lars S. Madsen, Vladyslav C. Usenko, Mikael Lassen, Radim Filip, and Ulrik L. Andersen.  Continuous variable quantum key distribution with modulated entangled states.  Nature Communications, 3(1):1083, Sep 2012.
  * [11] Tobias Gehring, Vitus Händchen, Jörg Duhme, Fabian Furrer, Torsten Franz, Christoph Pacher, Reinhard F. Werner, and Roman Schnabel.  Implementation of continuous-variable quantum key distribution with composable and one-sided-device-independent security against coherent attacks.  Nature Communications, 6(1):8795, Oct 2015.
  * [12] Han-Sen Zhong, Hui Wang, Yu-Hao Deng, Ming-Cheng Chen, Li-Chao Peng, Yi-Han Luo, Jian Qin, Dian Wu, Xing Ding, Yi Hu, Peng Hu, Xiao-Yan Yang, Wei-Jun Zhang, Hao Li, Yuxuan Li, Xiao Jiang, Lin Gan, Guangwen Yang, Lixing You, Zhen Wang, Li Li, Nai-Le Liu, Chao-Yang Lu, and Jian-Wei Pan.  Quantum computational advantage using photons.  Science, 370(6523):1460–1463, 2020.
  * [13] Lars S. Madsen, Fabian Laudenbach, Mohsen Falamarzi. Askarani, Fabien Rortais, Trevor Vincent, Jacob F. F. Bulmer, Filippo M. Miatto, Leonhard Neuhaus, Lukas G. Helt, Matthew J. Collins, Adriana E. Lita, Thomas Gerrits, Sae Woo Nam, Varun D. Vaidya, Matteo Menotti, Ish Dhand, Zachary Vernon, Nicolás Quesada, and Jonathan Lavoie.  Quantum computational advantage with a programmable photonic processor.  Nature, 606(7912):75–81, Jun 2022.
  * [14] Nicolas C. Menicucci, Peter van Loock, Mile Gu, Christian Weedbrook, Timothy C. Ralph, and Michael A. Nielsen.  Universal quantum computation with continuous-variable cluster states.  Phys. Rev. Lett., 97:110501, Sep 2006.
  * [15] Jun-ichi Yoshikawa, Shota Yokoyama, Toshiyuki Kaji, Chanond Sornphiphatphong, Yu Shiozawa, Kenzo Makino, and Akira Furusawa.  Invited Article: Generation of one-million-mode continuous-variable cluster state by unlimited time-domain multiplexing.  APL Photonics, 1(6):060801, 09 2016.
  * [16] Moran Chen, Nicolas C. Menicucci, and Olivier Pfister.  Experimental realization of multipartite entanglement of 60 modes of a quantum optical frequency comb.  Phys. Rev. Lett., 112:120505, Mar 2014.
  * [17] Warit Asavanant, Yu Shiozawa, Shota Yokoyama, Baramee Charoensombutamon, Hiroki Emura, Rafael N. Alexander, Shuntaro Takeda, Jun ichi Yoshikawa, Nicolas C. Menicucci, Hidehiro Yonezawa, and Akira Furusawa.  Generation of time-domain-multiplexed two-dimensional cluster state.  Science, 366(6463):373–376, 2019.
  * [18] Nathan Killoran, Thomas R. Bromley, Juan Miguel Arrazola, Maria Schuld, Nicolás Quesada, and Seth Lloyd.  Continuous-variable quantum neural networks.  Phys. Rev. Res., 1:033063, Oct 2019.
  * [19] Pere Mujal, Rodrigo Martínez-Peña, Johannes Nokkala, Jorge García-Beni, Gian Luca Giorgi, Miguel C. Soriano, and Roberta Zambrini.  Opportunities in quantum reservoir computing and extreme learning machines.  Advanced Quantum Technologies, 4(8):2100027, 2021.
  * [20] Herbert Jaeger.  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, 148(34):13, 2001.
  * [21] Wolfgang Maass, Thomas Natschläger, and Henry Markram.  Real-time computing without stable states: A new framework for neural computation based on perturbations.  Neural computation, 14(11):2531–2560, 2002.
  * [22] Mantas Lukoševičius and Herbert Jaeger.  Reservoir computing approaches to recurrent neural network training.  Computer Science Review, 3(3):127–149, 2009.
  * [23] F. Wyffels and B. Schrauwen.  A comparative study of reservoir computing strategies for monthly time series prediction.  Neurocomputing, 73(10):1958–1964, 2010.  Subspace Learning / Selected papers from the European Symposium on Time Series Prediction.
  * [24] Xiaowei Lin, Zehong Yang, and Yixu Song.  Short-term stock price prediction based on echo state networks.  Expert Systems with Applications, 36(3, Part 2):7313–7317, 2009\.
  * [25] Iulian Ilies, Herbert Jaeger, Olegas Kosuchinas, and Monserrat Rincon.  Stepping forward through echoes of the past: forecasting with echo state networks.  Technical report, 2007.
  * [26] Kohei Nakajima and Ingo Fischer.  Reservoir Computing: Theory, Physical Implementations, and Applications.  Springer, 2021.
  * [27] Fabian Triefenbach, Azarakhsh Jalalvand, Benjamin Schrauwen, and Jean-pierre Martens.  Phoneme recognition with large hierarchical reservoirs.  In J. Lafferty, C. Williams, J. Shawe-Taylor, R. Zemel, and A. Culotta, editors, Advances in Neural Information Processing Systems, volume 23. Curran Associates, Inc., 2010.
  * [28] Lin Wang, Zhigang Wang, and Shan Liu.  An effective multivariate time series classification approach using echo state network and adaptive differential evolution algorithm.  Expert Systems with Applications, 43:237–249, 2016.
  * [29] Gouhei Tanaka, Toshiyuki Yamane, Jean Benoit Héroux, Ryosho Nakane, Naoki Kanazawa, Seiji Takeda, Hidetoshi Numata, Daiju Nakano, and Akira Hirose.  Recent advances in physical reservoir computing: A review.  Neural Networks, 115:100–123, 2019.
  * [30] Daniel Brunner, Miguel C Soriano, Claudio R Mirasso, and Ingo Fischer.  Parallel photonic information processing at gigabyte per second data rates using transient states.  Nature communications, 4(1):1364, 2013.
  * [31] Kristof Vandoorne, Pauline Mechet, Thomas Van Vaerenbergh, Martin Fiers, Geert Morthier, David Verstraeten, Benjamin Schrauwen, Joni Dambre, and Peter Bienstman.  Experimental demonstration of reservoir computing on a silicon photonics chip.  Nature communications, 5(1):3541, 2014.
  * [32] Laurent Larger, Antonio Baylón-Fuentes, Romain Martinenghi, Vladimir S Udaltsov, Yanne K Chembo, and Maxime Jacquot.  High-speed photonic reservoir computing using a time-delay-based architecture: Million words per second classification.  Physical Review X, 7(1):011015, 2017.
  * [33] Guy Van Der Sande, Daniel Brunner, and Miguel C. Soriano.  Advances in photonic reservoir computing.  Nanophotonics, 6(3):561–576, 2017.
  * [34] Keisuke Fujii and Kohei Nakajima.  Harnessing disordered-ensemble quantum dynamics for machine learning.  Phys. Rev. Appl., 8:024030, Aug 2017.
  * [35] R. Martínez-Peña, J. Nokkala, G. L. Giorgi, R. Zambrini, and M. C. Soriano.  Information processing capacity of spin-based quantum reservoir computing systems.  Cognitive Computation, Oct 2020.
  * [36] Johannes Nokkala, Rodrigo Martínez-Peña, Gian Luca Giorgi, Valentina Parigi, Miguel C. Soriano, and Roberta Zambrini.  Gaussian states of continuous-variable quantum systems provide universal and versatile reservoir computing.  Communications Physics, 4(1):53, Mar 2021.
  * [37] Michele Spagnolo, Joshua Morris, Simone Piacentini, Michael Antesberger, Francesco Massa, Andrea Crespi, Francesco Ceccarelli, Roberto Osellame, and Philip Walther.  Experimental photonic quantum memristor.  Nature Photonics, 16(4):318–323, 2022.
  * [38] Jiayin Chen, Hendra I. Nurdin, and Naoki Yamamoto.  Temporal information processing on noisy quantum computers.  Phys. Rev. Appl., 14:024065, Aug 2020.
  * [39] Antonio Sannia, Rodrigo Martínez-Peña, Miguel C Soriano, Gian Luca Giorgi, and Roberta Zambrini.  Dissipation as a resource for quantum reservoir computing.  arXiv preprint arXiv:2212.12078, 2022.
  * [40] Sanjib Ghosh, Andrzej Opala, Michał Matuszewski, Tomasz Paterek, and Timothy C. H. Liew.  Quantum reservoir processing.  npj Quantum Information, 5(1):35, Apr 2019.
  * [41] Guillem Llodrà, Christos Charalambous, Gian Luca Giorgi, and Roberta Zambrini.  Benchmarking the role of particle statistics in quantum reservoir computing.  Advanced Quantum Technologies, 6(1):2200100, 2023.
  * [42] Rodrigo Martínez-Peña, Gian Luca Giorgi, Johannes Nokkala, Miguel C. Soriano, and Roberta Zambrini.  Dynamical phase transitions in quantum reservoir computing.  Phys. Rev. Lett., 127:100502, Aug 2021.
  * [43] Johannes Nokkala, Rodrigo Martínez-Peña, Roberta Zambrini, and Miguel C. Soriano.  High-performance reservoir computing with fluctuations in linear networks.  IEEE Transactions on Neural Networks and Learning Systems, 33(6):2664–2675, 2022.
  * [44] Pere Mujal, Rodrigo Martínez-Peña, Gian Luca Giorgi, Miguel C. Soriano, and Roberta Zambrini.  Time-series quantum reservoir computing with weak and projective measurements.  npj Quantum Information, 9(1):16, Feb 2023.
  * [45] Jorge García-Beni, Gian Luca Giorgi, Miguel C. Soriano, and Roberta Zambrini.  Scalable photonic platform for real-time quantum reservoir computing.  Phys. Rev. Appl., 20:014051, Jul 2023.
  * [46] Fangjun Hu, Gerasimos Angelatos, Saeed A Khan, Marti Vives, Esin Türeci, Leon Bello, Graham E Rowlands, Guilhem J Ribeill, and Hakan E Türeci.  Tackling sampling noise in physical systems for machine learning applications: Fundamental limits and eigentasks.  arXiv preprint arXiv:2307.16083, 2023.
  * [47] W. D. Kalfus, G. J. Ribeill, G. E. Rowlands, H. K. Krovi, T. A. Ohki, and L. C. G. Govia.  Hilbert space as a computational resource in reservoir computing.  Phys. Rev. Res., 4:033007, Jul 2022.
  * [48] M. C. Soriano, S. Ortín, D. Brunner, L. Larger, C. R. Mirasso, I. Fischer, and L. Pesquera.  Optoelectronic reservoir computing: tackling noise-induced performance degradation.  Opt. Express, 21(1):12–20, Jan 2013.
  * [49] M.C. Soriano, Silvia Ortin, Lars Keuninckx, Lennert Appeltant, Jan Danckaert, L. Pesquera, and Guy Van Der Sande.  Delay-based reservoir computing: Noise effects in a combined analog and digital implementation.  IEEE Transactions on Neural Networks and Learning Systems, 26(2):388–393, February 2015.
  * [50] Pere Mujal, Rodrigo Martínez-Peña, Gian Luca Giorgi, Miguel C. Soriano, and Roberta Zambrini.  Time-series quantum reservoir computing with weak and projective measurements.  npj Quantum Information, 9(1):16, Feb 2023.
  * [51] Kohei Nakajima and Ingo Fischer.  Reservoir Computing: Theory, Physical Implementations, and Applications.  Springer Singapore, 2021.
  * [52] Lyudmila Grigoryeva and Juan-Pablo Ortega.  Echo state networks are universal.  Neural Networks, 108:495–508, 2018.
  * [53] R. Medeiros de Araújo, J. Roslund, Y. Cai, G. Ferrini, C. Fabre, and N. Treps.  Full characterization of a highly multimode entangled state embedded in an optical frequency comb using pulse shaping.  Phys. Rev. A, 89:053828, May 2014.
  * [54] Jonathan Roslund, Renné Medeiros de Araújo, Shifeng Jiang, Claude Fabre, and Nicolas Treps.  Wavelength-multiplexed quantum networks with ultrafast frequency combs.  Nature Photonics, 8(2):109–112, Feb 2014.
  * [55] Y. Cai, J. Roslund, G. Ferrini, F. Arzani, X. Xu, C. Fabre, and N. Treps.  Multimode entanglement in reconfigurable graph states using optical frequency combs.  Nature Communications, 8(1):15645, Jun 2017.
  * [56] Tiphaine Kouadou, Francesca Sansavini, Matthieu Ansquer, Johan Henaff, Nicolas Treps, and Valentina Parigi.  Spectrally shaped and pulse-by-pulse multiplexed multimode squeezed states of light.  APL Photonics, 8(8), 2023.
  * [57] Lorenz Butschek, Akram Akrout, Evangelia Dimitriadou, Alessandro Lupo, Marc Haelterman, and Serge Massar.  Photonic reservoir computer based on frequency multiplexing.  Optics Letters, 47(4):782–785, 2022.
  * [58] J Nokkala, F Arzani, F Galve, R Zambrini, S Maniscalco, J Piilo, N Treps, and V Parigi.  Reconfigurable optical implementation of quantum complex networks.  New Journal of Physics, 20(5):053024, may 2018.
  * [59] Albert Cabot, Fernando Galve, Víctor M. Eguíluz, Konstantin Klemm, Sabrina Maniscalco, and Roberta Zambrini.  Unveiling noiseless clusters in complex quantum networks.  npj Quantum Information, 4(1):57, Nov 2018.
  * [60] Paul Renault, Johannes Nokkala, Ganaël Roeland, Nicolas Joly, Roberta Zambrini, Sabrina Maniscalco, Jyrki Piilo, Nicolas Treps, and Valentina Parigi.  Experimental optical simulator of reconfigurable and complex quantum environment, 2023.
  * [61] S. Ortín, M. C. Soriano, L. Pesquera, D. Brunner, D. San-Martín, I. Fischer, C. R. Mirasso, and J. M. Gutiérrez.  A unified framework for reservoir computing and extreme learning machines based on a single time-delayed neuron.  Scientific Reports, 5(1):14945, Oct 2015.
  * [62] J. Doyne Farmer and John J. Sidorowich.  Predicting chaotic time series.  Phys. Rev. Lett., 59:845–848, Aug 1987.
  * [63] Herbert Jaeger and Harald Haas.  Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication.  Science, 304(5667):78–80, 2004.
  * [64] Tomoyuki Kubota, Hirokazu Takahashi, and Kohei Nakajima.  Unifying framework for information processing in stochastically driven dynamical systems.  Phys. Rev. Res., 3:043135, Nov 2021.
  * [65] Joni Dambre, David Verstraeten, Benjamin Schrauwen, and Serge Massar.  Information processing capacity of dynamical systems.  Scientific reports, 2:514, 2012.
  * [66] Masanobu Inubushi and Kazuyuki Yoshimura.  Reservoir computing beyond memory-nonlinearity trade-off.  Scientific reports, 7:10199, 2017.
  * [67] Michael C. Mackey and Leon Glass.  Oscillation and chaos in physiological control systems.  Science, 197(4300):287–289, 1977.
  * [68] Daniel Brunner, Miguel C Soriano, and Guy Van der Sande.  Photonic reservoir computing.  De Gruyter, 2019.
  * [69] Michele Spagnolo, Joshua Morris, Simone Piacentini, Michael Antesberger, Francesco Massa, Andrea Crespi, Francesco Ceccarelli, Roberto Osellame, and Philip Walther.  Experimental photonic quantum memristor.  Nature Photonics, 16(4):318–323, 2022.
  * [70] Samuel L. Braunstein.  Squeezing as an irreducible resource.  Phys. Rev. A, 71:055801, May 2005.
  * [71] Gianfranco Cariolaro and Gianfranco Pierobon.  Reexamination of Bloch-Messiah reduction.  Phys. Rev. A, 93:062115, Jun 2016.
  * [72] Roger A. Horn and Charles R. Johnson.  Matrix Analysis.  Cambridge University Press, 1985.
  * [73] Zoran Konkoli.  On reservoir computing: from mathematical foundations to unconventional applications.  In Advances in unconventional computing, pages 573–607. Springer, 2017.