Experimental memory control in continuous-variable optical quantum reservoir computing

The quantum resource

The core quantum resource is a multimode optical network generated via parametric down-conversion (PDC) in a nonlinear waveguide pumped by a pulsed laser, as shown as in Fig. 1. The interaction Hamiltonian describing the generation of signal and idler fields of wavelengths λ_S, λ_I reads

where J(λ_S, λ_I), the so-called joint spectral amplitude (JSA), can be written as^45,46,47

The term p(λ_S, λ_I) is the complex pump spectral shape, while ϕ(λ_S, λ_I) is the phase-matching function, determined by the characteristics of the nonlinear medium. The multimode nature of the generated quantum system can be retrieved by performing a Schmidt decomposition of the JSA J(λ_S, λ_I) = ∑_kr_kh_k(λ_S)g_k(λ_I), where h_k(λ_S) and g_k(λ_I) are two sets of orthonormal modes, spectral profiles of signal and idler channels, and r_k are called the Schmidt coefficients. In case of degenerate signal and idler fields, as in our experimental set-up, h_k = g_k; and the Hamiltonian in equation (1) can be written as \(\widehat{H}\propto {{\rm{\sum }}}_{k}{r}_{k}\,{\widehat{A}}_{k}^{\dagger 2}+{\rm{h}}.{\rm{c}}.\), with \({\widehat{A}}_{k}^{\dagger }\) broadband creation operators associated to the modes h_k. Acting on the vacuum, this Hamiltonian generates a tensor product of single-mode squeezed vacuum states, one in each mode h_k(λ_S), also called supermodes. The used set-up can generate up to 40 Schmidt modes approximated by Hermite–Gauss polynomials as shown in Fig. 1.

The generated quantum state is Gaussian and fully characterized by the quadrature covariance matrix σ, measured via homodyne detection. Mode-selective homodyne detection enables probing in arbitrary mode bases. In the supermode basis, σ is diagonal; in any other basis ξ_j(λ_S) related by a unitary U, it transforms as \({{\boldsymbol{\sigma }}}^{{\prime} }={{\boldsymbol{S}}}_{{\boldsymbol{U}}}^{T}{\boldsymbol{\sigma }}{{\boldsymbol{S}}}_{{\boldsymbol{U}}}\), where S_U is the symplectic transformation induced by U (Supplementary Section B). The matrix \({{\boldsymbol{\sigma }}}^{{\prime} }\) is non-diagonal and the off-diagonal terms are signatures of entanglement correlations between the measured modes. Measuring on such basis allows to experimentally access entanglement correlations for quantum protocols^35,48. The resource can then be pictured via an entangled network as shown in Fig. 2.

The phase-matching function term of the JSA in equation (2) is fixed and determined by the material and waveguide properties, and only the pump-envelope function can be modified in a reconfigurable way via mode-shaping techniques. This modifies the entanglement correlation between nodes as shown in Fig. 2. Therefore, we can explicitly write the dependence of the covariance matrix \({{\boldsymbol{\sigma }}}^{{\prime} }\) on the pump p and the chosen basis U as \({{\boldsymbol{\sigma }}}^{{\prime} }={{\bf{S}}}_{{\bf{U}}}^{T}\,{\boldsymbol{\sigma }}({\bf{H}}({\bf{p}}))\,{{\bf{S}}}_{{\bf{U}}}\).

In the following, we exploit the control on the pump-envelope function to encode information in the quantum system in the configuration of a fixed measurement basis U. The expectation value of general observables \({O}_{m}=\langle \Delta {\widehat{x}}_{i}\Delta {\widehat{x}}_{j}\rangle\) that will be used in the reservoir task can be retrieved from the covariance matrix \({{\boldsymbol{\sigma }}}^{{\prime} }(U,p)\). Given the pump shape, we have accurate numerical models of the process from the nonlinear interaction to the observables^45,46,47, setting a Digital Twin of the experimental set-up as fully explained in Supplementary Section A.

Information encoding and memory control of the quantum reservoir

Tuning the pump phase profile, the pump-envelope function can be written as

where the f_i(λ_P) are a set of N pump modes and the \({\boldsymbol{\delta }}={({\delta }_{i})}_{i=1}^{N}\in {{\mathbb{R}}}^{N}\) are tunable phase parameters that allow precise shaping of the JSA to easily encode information. The pump input modes f_i that we consider are N frequency windows, delimiting a portion of our Gaussian pump profile f(λ_P), that we can precisely control with a spatial light modulator, for example. In the simple case of N = 1, when a unique global phase parameter is considered, an EOM is used as shown in Fig. 1. To perform the measurement, a set of n modes ξ_j(λ_S) is selected by the basis choice U via the local oscillator shaping in homodyne. The tailoring of N pump spectral components and the choice of the measurement basis composed by n elements are shown in the upper part of Fig. 2 (more details in Supplementary Section B).

Our system is endowed with memory through a real-time feedback mechanism: the results of previous measurements are fed back into the system to influence future states. Specifically, if we consider a discrete-time sequence of inputs {s_k}, where k denotes the timestep, the pump parameters are updated according to

Here i ∈ 1, …, N, α_i and β_i are real-valued coefficients, and M_im is a feedback mask (a real matrix) that determines how previous measurement outcomes influence the current pump parameters. This determines the fading memory of the QRC. To summarize, N is the number of different input components that we encode in a single instance by manipulating several components of the pump spectrum: N > 1 can encode multidimensional inputs as well as enrich the encoding of a scalar input, repeated across the pump spectrum with different parameters (α_i, β_i, M_im). The quantities \({O}_{m}^{(k-1)}={\langle \Delta {\widehat{x}}_{i}\Delta {\widehat{x}}_{j}\rangle }^{(k-1)}\) are the values of the mth observable measured at timestep k − 1, and are derived from a subset of elements of the covariance matrix (i, j ∈ 1, …, 2n). These are measured by accumulating data on the homodyne signal. The number n of measured modes determines then the size of the output layer of our reservoir. The reservoir computing is described by

The reservoir is trained by adjusting the matrix w and parameter b to minimize the error in the chosen learning task. The scheme of the quantum reservoir with feedback is shown in the top part of Fig. 1. While we here focus on the case of feedback restricted to the last output at each cycle, the fading memory can be extended re-injecting a longer combination of past signals.

QRC tasks

We evaluate the performance of our system across a range of benchmark tasks considering two input encoding: the simple, experimentally implemented scheme with N = 1, and the more general encoding with N > 1. We begin with the XOR task, which assesses the system’s nonlinearity and short-term memory while requiring minimal resources. We then consider the memory task to characterize the system’s fading memory property. Finally, we turn to more complex benchmarks, including the double-scroll and parity check tasks, which challenge the system’s scalability and capacity to handle richer, higher-order nonlinear dynamics.

Experimental control of memory via global phase

The phase encoding via the EOM, as shown in Fig. 1, corresponds to N = 1, and the voltage applied changes linearly the pump phase δ:

where f(λ_P) is the natural Gaussian profile of the pump. As explained in Methods, with this type of encoding, the expressivity of the observables does not increase with the number of measured modes n; therefore, we use n = 1 as well. This means that we measure only the first supermode h₁ and we can use as observables {O₁, O₂, O₃} the set \(\{\langle \Delta {\widehat{q}}^{2}\rangle ,\langle \Delta {\widehat{p}}^{2}\rangle ,\langle \Delta \widehat{q}\Delta \widehat{p}\rangle \sim \langle \Delta {(\frac{\widehat{q}+\widehat{p}}{\sqrt{2}})}^{2}\rangle \}\). These three observables are sufficient for implementing several tasks, as they express three different nonlinear dependencies with respect to δ as shown in the lower part of Fig. 1. The first two panels show the three behaviours experimentally, in the case of a low-noise realization and a more noisy one, which is representative of the average noise in measurements. Assuming a Gaussian noise model, the typical variance of observables is recovered via a least square method. The result is used in the Digital Twin to simulate the expected behaviour as shown in the last two panels.

The unique non-zero term M that we use for the feedback mask is the one corresponding to \({\langle \Delta {\widehat{q}}^{2}\rangle }^{(k-1)}\). At each timestep k, we set the global phase depending on the input s_k and the feedback as

M, α and β are randomly chosen within ranges that ensure both the input and the feedback contribute comparably to the phase and such to explore a nontrivial region of the observable space.

We begin with a binary task, the temporal XOR, where the target output is \({\widehat{y}}_{k}={s}_{k}\bigoplus {s}_{k-1}\) and the input is a random sequence of values s_k ∈ {0, 1}. This task requires both nonlinearity and memory of the previous encoded input. The reservoir is trained on a set of inputs by applying linear regression to the measured observables to maximize accuracy (Supplementary Section C).

The performance of the reservoir is then evaluated using the accuracy again, but on a test sequence of new data. In Fig. 3a, we present an example of the test target values and the corresponding experimental prediction. In Fig. 3b,c, we show experimental and simulated accuracies as a function of training set size, for different levels of Gaussian noise. Experimental instabilities can introduce Gaussian noise in the measurement of the observables (Fig. 1), which can be mitigated by averaging over repeated measurements. These instabilities arise from mechanical and thermal fluctuations and are not intrinsic limitations of the protocol, but rather reflect the current performance of the experimental set-up. We see that the training saturates for set sizes around 100 and that noise affects the performance. Nevertheless, the reservoir could achieve experimentally test accuracies of 98 ± 1% for realistic levels of noise, with only 70 training steps. The Digital Twin behaves similarly to the experiment, as shown in Fig. 3c.

A way to further test the memory of the system is through the linear memory task, in which the reservoir is trained to reproduce past entries of the input series. The target function at timestep k is the input encoded τ steps in the past, \({\widehat{y}}_{k}={s}_{k-\tau }\), the input s_k ∈ [−1, 1] being a random uniformly distributed sequence. Here the three previous observables are not sufficient. Therefore, we resort to spatial multiplexing, which is equivalent to use R different reservoirs at the same time. This is a classical reservoir computing technique used to enhance the expressivity of the system (see Methods for details), whereas the general encoding introduced in the next section leverages spectral multiplexing, an intrinsic resource of the quantum system. In Fig. 3d, we present an example of target and experimentally predicted values. In Fig. 3e, instead, we show the agreement between the experimental and simulated capacities (see Supplementary Section C for the definition) on test sets when varying τ and the number R of reservoirs. The simulations agree with the experimental results, confirming the advantage of multiplexing (better capacities for higher R) and the fading memory of the system (capacity gradually decreasing for higher delays).

Given the optimal agreement with the experiment, we use the Digital Twin to test the protocol on a more complex task. Concretely, we perform forecasting of chaotic time series such as the double-scroll electronic circuit⁴⁹ and also (Supplementary Section H) the Lorenz system. Since the input vectors are three-dimensional, for the global phase encoding, we require at least three spatially multiplexed reservoirs (one per input signal) with additional reservoirs possible to enhance expressivity. The target function is the input one timestep in the future: \({\widehat{{\bf{y}}}}_{k}={{\bf{s}}}_{k+1}\). Concerning the double-scroll, in Fig. 4a, we present an example of target and predicted values, obtained with simulations with realistic noise values. In Supplementary Fig. 2, we present three plots (corresponding to V₁, V₂ and I) where the capacity is calculated while varying the training size and the number of reservoirs R. We prove that even with the experimental noise, we saturate the prediction capacity for relatively small values of R (below 15) and of training size (below 400).

Scaling the expressivity and boosting the memory with general encoding

In the global phase encoding, the expressivity is limited for any number of measured modes n > 1 (Methods). Therefore, to take advantage of the multimode nature of the system, we diversify the phase profile of the pump (N > 1). This means exploiting fully the phase encoding from equation (4), as shown in Fig. 4. In Fig. 4b, with this general encoding (N = 5), we achieve perfect capacities (100%) for a parity check task for τ ∈ {1, 2, 3} and prove the boost in memory when n, the number of measurement modes, is increased. Furthermore, as shown in Fig. 4c, the expressivity (see Supplementary Section C for definition) increases with n for N > 1. Specifically, for N > 5, we saturate a polynomial scaling with n as reported in ref. ²⁵: the expressivity scales with n(n + 1)/2, which is also the dimension of the observables matrix (considering only \(\widehat{q}\) observables). This means that with a sufficiently large N, all the terms in our covariance matrix are independent. We highlight the increased diversity in behaviours of the observables in Fig. 4d. We also exploit the multidimensionality of the input (it can be a vector of size up to N), to implement again the double-scroll task, without multiplexing. The performances of the global phase encoding with multiplexing (R = 15 reservoirs) and the general encoding with N, n = 9, 9 are similar, but the gain in scalability is clear as with the general encoding use only one system, instead of multiplexing in 15 reservoirs.