# Fully analog photonic deep reservoir computer based on frequency
multiplexing

Alessandro<EMAIL_ADDRESS>Laboratoire d’Information Quantique,
CP 224, Universitè libre de Bruxelles, Av. F. D. Roosevelt 50, B-1050,
Bruxelles, Belgium  Enrico Picco  Laboratoire d’Information Quantique, CP 224,
Universitè libre de Bruxelles, Av. F. D. Roosevelt 50, B-1050, Bruxelles,
Belgium  Marina Zajnulina  Laboratoire d’Information Quantique, CP 224,
Universitè libre de Bruxelles, Av. F. D. Roosevelt 50, B-1050, Bruxelles,
Belgium  Serge Massar  Laboratoire d’Information Quantique, CP 224, Universitè
libre de Bruxelles, Av. F. D. Roosevelt 50, B-1050, Bruxelles, Belgium

###### Abstract

Reservoir computers (RC) are randomized recurrent neural networks well adapted
to process time series, performing tasks such as nonlinear distortion
compensation or prediction of chaotic dynamics. Deep reservoir computers
(deep-RC), in which the output of one reservoir is used as input of the next
reservoir, can lead to improved performance because, as in other deep
artificial neural networks, the successive layers represent the data in more
and more abstract ways. We present a fiber-based photonic implementation of a
two-layer deep-RC based on frequency multiplexing. The two RC layers are
encoded in two frequency combs propagating in the same experimental setup. The
connection between layers is fully analog and does not require any digital
processing. We find that the deep-RC outperforms traditional RC by up to two
orders of magnitude on two benchmark tasks. This work thus paves the way
towards using fully analog photonic neuromorphic computing for complex
processing of time series, while avoiding costly analog-to-digital and
digital-to-analog conversions.

##  1 Introduction

Artificial Intelligence is probably the most disruptive new technology to
emerge during the first decades of the XXIst century. Its success is based on
the use of deep neural networks in which multiple layers of artificial neurons
are connected in a feed forward architecture [4, 19]. Recent advances include,
for instance, image classification and analysis [28], game playing [31],
protein structure prediction[3, 16], chat bots that simulate human
conversation such as ChatGPT and Bing[26, 5], and more.

Artificial neural networks are fundamentally analog systems simulated on a
digital computer. Thus, it seems very attractive to replace the digital
simulation by analog hardware, as this could result in considerable energy
savings. Photonics is particularly attractive for such analog neural networks
because of its potential for very high speed (see e.g. [39, 8]), parallelism
(see e.g. [22, 29]), possibility of implementing spiking networks (see e.g.
[7, 15]), and low energy consumption per operation (see e.g. [11]). The
importance of deep neural networks for complex applications has lead to
several demonstrations of deep photonic networks [30, 21, 40, 2].

It is also possible to use physical substrates that process information in
ways very different from traditional artificial neural networks, such as
extreme learning machines[13] and reservoir computers[14], see e.g. the
reviews [33, 37, 23]. These physical substrates for computing are very
attractive because of the augmented freedom in the design that they offer, for
instance, in the combination of several systems to constitute deep networks.
However, the training and implementation of deep physical networks provide new
challenges which are only starting to be addressed [38, 25].

Reservoir computers (RC), which are the topic of this work, are randomized
recurrent neural networks (RNN) in which the recurrency is provided by a
(simulated or physical) high dimensional nonlinear dynamical system called
“reservoir” [14]. RCs have been successfully implemented in analog systems
including photonics, electronics, spintronics, mechanics, biology and more
(see e.g. [33]). Many of the photonic implementations of RC used a delay loop
and a single dynamical node, encoding multiple neurons using time multiplexing
[1]. Even though high speed implementations have been realised using this
approach [18], the time multiplexing represents an inherent slowdown.
Alternative approaches include spatial multiplexing, e.g. using free space
optics [27], multimode fibers [32], as well as hybrid temporal/spatial
approaches [24].

As for other neural network architectures, assembling several RCs in a deep
architecture enhances performance. Deep RCs were first used in [34] and
studied in more depth in [10], where it is shown that the serial connection
among different RC layers enhances the system performance by enriching its
dynamics. Different ways of combining photonic reservoirs into networks are
compared in [9]. Motivated by these works, a first experimental implementation
of deep-RC is reported in [25], showing significant improvement in performance
with increasing the number of layers. However, in this work each reservoir was
implemented using the time multiplexing architecture, which is not optimal in
terms of computing speed, and, more importantly, the connection between
reservoirs was implemented digitally. The latter is also the case in the
related work [38]. Ref. [20] proposes an architecture for a deep reservoir
based on time delay architecture with analog connection between layers.

Here we report a deep reservoir computer in which the connection between
layers is fully analog and does not require storage or processing on a digital
computer. Our experiment is based on a recently reported reservoir computer in
which the neuron signals are encoded in the amplitudes of a frequency comb,
and mixing between neurons is realised by electro-optic phase modulators [6].
This architecture allows for a relatively easy to realize optical output
layer, as weights can be applied on comb lines by a programmable spectral
filter, and the nonlinear summation of the weighted neuron signals can be
executed by a photodiode, which measures the total intensity of the weighted
frequency comb and introduces a quadratic nonlinearity. This technique,
already employed in [6] to generate the output signals with optical weighting,
in the present work allows us to use the output of a reservoir as input to a
second one without leaving the analog domain. In the present work we also
fully exploit the frequency degree of freedom of the light by using the same
hardware for implementing multiple reservoirs simultaneously, each one working
in a different frequency band. In particular, here we report two simultaneous
RC computations and we demonstrate that combining the two computation in a
deep fashion improves performance compared to using the two reservoirs in
parallel without interconnections. We test two strategies for optimizing the
interconnections between two layers in the deep configuration. In the first,
simpler, approach we only adjust the strength of the connection; while in the
second approach we optimise the connections using the Covariance Matrix
Adaptation Evolution Strategy (CMA-ES)[12]. Surprisingly we find that both
approaches yield comparable results.

In Sec. 2 we present the algorithms, the experimental setup and the
benchmarking methods; in Sec. 3 we present and discuss results; finally, in
Sec. 4 we present conclusions and outlooks for this work.

##  2 Methods

###  2.1 Algorithms

Figure 1: Left panel: standard reservoir computing scheme. Right panel: deep
reservoir computing scheme. The weights in black are fixed, while the weights
in red are trained.

####  2.1.1 Reservoir computing

A reservoir computer (RC, see left panel of Fig.1) [14] is a recurrent neural
network composed of three layers: the input layer, the reservoir layer and the
output layer. Only the output weights are trained, while the input and
internal weights are fixed and not trained.

The experimental system is based on the frequency multiplexing RC scheme
described in [6]. Neurons are encoded in the complex amplitudes of the lines
of a frequency comb and neuronal interconnections are realized via frequency-
domain interference that makes comb line exchange power. The electric field at
any point in the reservoir can thus be expressed as

$E(t)=\sum_{k}x_{k}(t)\exp\left(i(\omega+k\Omega)t\right),$ (1)

where $\omega$ is the center frequency of the comb, $\Omega$ the frequency
spacing between comb lines, and $x_{k}(t)$ are the slowly varying amplitudes
of the comb lines which encode neuron information. To describe more
conveniently the RC application, we restrict to the $N$ most central lines of
the comb, which are the ones encoding information. Moreover, we group the
amplitude of these lines in a $N$ dimensional complex vector $\mathbf{x}_{n}$
that evolves in slow, discrete, time $n$. The discrete timescale corresponds
to the discrete evolution of the RC states.

The RC based on frequency multiplexing uses nonlinear input and output layers,
and a linear reservoir (which is a powerful architecture, as demonstrated in
[35]). It can be described by the evolution equations:

$\displaystyle\mathbf{x}_{n}$
$\displaystyle=\mathbf{W}_{\text{}}\cdot\mathbf{x}_{n-1}+\mathbf{W}_{\text{in}}\cdot
f_{\text{in}}\left(u_{n}\right),$ (2) $\displaystyle{y}_{n}$
$\displaystyle=\left|\mathbf{W}_{\text{out}}^{+}\cdot\mathbf{x}_{n}\right|^{2}-\left|\mathbf{W}_{\text{out}}^{-}\cdot\mathbf{x}_{n}\right|^{2},$
(3)

where $u_{n}$ (a real scalar) is the input signal to the reservoir at timestep
${n}$, and ${y}_{n}$ (a real scalar) is the output signal of the reservoir at
timestep $n$, $\mathbf{W}_{\text{}}$ is a complex $N\times N$ matrix
representing the internal connections of the reservoir,
$\mathbf{W}_{\text{in}}$ is a complex $N$ dimensional vector representing the
input-to-reservoir connections, and $\mathbf{W}_{\text{out}}^{+}$ and
$\mathbf{W}_{\text{out}}^{-}$ are $N\times N$ diagonal matrices with positive
real coefficients representing the output weights. In this notation,
$\mathbf{W}_{\text{out}}^{-}\cdot\mathbf{x}_{n}$ and
$\mathbf{W}_{\text{out}}^{+}\cdot\mathbf{x}_{n}$ are complex $N$ dimensional
vectors representing the the comb line amplitudes after the application of
negative and positive weights respectively.

The input signal is supplied through a Mach-Zehnder modulator operating in the
negative quadrature point, hence the input nonlinearity $f_{\text{in}}$ is
given by the modulator transfer function:

$f_{\text{in}}(u)=E_{0}\cdot\sin\left(\gamma\cdot u\right),$ (4)

where $E_{0}$ represents the input radiation amplitude and $\gamma$ is the
driving strength of the electrical signal to the modulator.

The vector norm square $|\cdot|^{2}$ in Eq. (3) represents the output
nonlinearity given by the quadratic response of photodiodes. Note that the
output ${y}_{n}$ can also be computed by first measuring the intensity of each
comb line independently, and then performing the multiplication by
$(\mathbf{W}_{\text{out}}^{+})^{2}$ or $-(\mathbf{W}_{\text{out}}^{-})^{2}$
according to whether the weight is positive or negative. The weights are
optimised using ridge regression so that the output $y_{n}$ approximates as
well as possible the desired output.

####  2.1.2 Deep Reservoir Computing

A deep reservoir computer (deep-RC, see right panel of Fig. 1) is a collection
of RC layers connected in series. The deep-RC output signal is a linear
combination of neuron values of each reservoir. The hierarchy introduced by
the serial connection enhances the network performance because the different
reservoirs can have independent dynamics, thus enriching the states of the
full deep-RC.

The deep-RC composed of $N_{\text{layers}}$ layers, each one comprising $N$
neurons, as implemented in our system, is described by the set of equations:

$\displaystyle\mathbf{x}_{n}^{\text{(1)}}$
$\displaystyle=\mathbf{W}_{\text{}}^{\text{(1)}}\cdot\mathbf{x}_{{n}-1}^{\text{(1)}}+\mathbf{W}_{\text{in}}^{\text{(1)}}\cdot
f_{\text{in}}\left(u_{n}\right),$ (5)
$\displaystyle\mathbf{x}_{n}^{\text{(i)}}$
$\displaystyle=\mathbf{W}_{\text{}}^{\text{(i)}}\cdot\mathbf{x}_{{n}-1}^{\text{(i)}}+\mathbf{W}_{\text{in}}^{\text{(i)}}\cdot
f_{\text{in}}\left(u_{n}^{(i)}\right),\quad$ $\displaystyle i$
$\displaystyle=2,\ldots,N_{\text{layers}}$ (6) $\displaystyle u_{n}^{(i+1)}$
$\displaystyle=\left|\mathbf{W}_{\text{out}}^{\text{(i)}}\cdot\mathbf{x}_{n}^{\text{(i)}}\right|^{2},\quad$
$\displaystyle i$ $\displaystyle=1,\ldots,N_{\text{layers}}-1$ (7)
$\displaystyle{y}_{n}^{(\text{T})}$
$\displaystyle=\left|{\mathbf{W}_{\text{out}}^{\text{(\text{T+})}}}\cdot\mathbf{x}_{n}^{\text{(\text{T})}}\right|^{2}-\left|{\mathbf{W}_{\text{out}}^{\text{(\text{T-})}}}\cdot\mathbf{x}_{n}^{\text{(\text{T})}}\right|^{2},$
(8)

where the superscript $(i)$, $1\leq i\leq N_{\text{layers}}$, identifies the
reservoir layer. Here, as before, $\mathbf{W}_{\text{}}^{\text{(i)}}$ is a
complex $N\times N$ matrix representing the internal connections of the $i$-th
reservoir layer, $\mathbf{W}_{\text{in}}^{\text{(i)}}$ is a complex $N$
dimensional vectors representing the input connections of the $i$-th reservoir
layer, $\mathbf{W}_{\text{out}}^{\text{(i)}}$ is a $N\times N$ diagonal
matrices with positive real coefficients representing the output connections
of the $i$-th layer. In our current photonic implementation,
$N_{\text{layers}}=2$, but the equations generalise easily to more layers. The
first reservoir layer is driven by the input time series $u_{n}$, while the
next reservoir layers are driven by the output $u_{n}^{(i)}$ of the previous
layer, see Eq. (7). Note that in our implementation the connections among
consecutive layers are intermediated by only positive weights, contained in
the diagonal of $\mathbf{W}_{\text{out}}^{\text{(i)}}$, which is why there is
only a single term on the right hand side of Eq. (7).

The output of the deep-RC, ${y}_{n}^{(\text{T})}$, is obtained by combing the
states from both layers, i.e. by multiplying the lines of both combs by output
weights, and measuring the sum with a photodiode. To this purpose, we have
defined
$\mathbf{x}_{n}^{\text{(\text{T})}}=\left(\mathbf{x}_{n}^{\text{(1)}},\mathbf{x}_{n}^{\text{(2)}},\ldots,\mathbf{x}_{n}^{(N_{\text{layers}})}\right)$,
as the complex vector of size $N_{\text{layers}}\cdot N$ representing the full
deep-RC state at timestep $n$, and
$\mathbf{W}_{\text{out}}^{\text{(\text{T+})}}$ and
$\mathbf{W}_{\text{out}}^{\text{(\text{T-})}}$ as the diagonal
$(N_{\text{layers}}\cdot N)\times(N_{\text{layers}}\cdot N)$ matrices with
positive real coefficients representing the deep-RC output weights. The output
weights are optimised using ridge regression.

Note that the interconnection between successive layers (say layers $i$ and
$i+1$) is determined by $3N$ real parameters: the $N$ positive real elements
of the diagonal matrix $\mathbf{W}_{\text{out}}^{\text{(i)}}$, and both the
real and imaginary parts of the $N$ elements of the vector
$\mathbf{W}_{\text{in}}^{\text{(i+1)}}$. Of these, only the $N$ elements of
the diagonal matrix $\mathbf{W}_{\text{out}}^{\text{(i)}}$ can be tuned in our
experimental setup. This is to be compared with the proposal of [10], in which
the interconnection is given by a $N\times N$ random matrix, whose spectral
radius is tuned.

###  2.2 Experimental setup

Figure 2: Experimental setup. Optical connections are in blue, electrical
connections in red. MZM: Lithium Niobate Mach-Zehnder modulator; AWG:
arbitrary waveform generator; C: fiber couplers; EDFA: Erbium-doped-fiber
amplifier; PM: phase modulator; RF source: radio frequency source at frequency
$\Omega$; RF AMP: radio frequency amplifiers; PSF: programmable spectral
filter; PD: photodiode; ES: electric switch. Figure 3: Normalized spectral
power of the radiation as measured at the output of the fiber loop, after
coupler C3. Red markers indicate the input wavelengths $\lambda_{1}=1550.2$ nm
and $\lambda_{2}=1555.4$ nm.

##### Experimental Setup.

The experimental system is based on what described in [6], modified in such a
way to support two RC computations at the same time.

Fig. 2 reports the schematic of the experiment. All fiber connections and
couplers are single-mode and polarization-maintaining. We employ two
continuous-wave laser sources (CW source 1 and CW source 2) at wavelengths
$\lambda_{1}=1550.2$ nm and $\lambda_{2}=1555.4$ nm. The two laser outputs are
modulated by two Mach-Zehnder modulators (MZM 1 and MZM 2). Both MZMs are
biased to operate in the negative quadrature point (bias controllers are not
shown in Fig. 2). The transfer functions of MZM 1 and MZM 2 define the input
nonlinearities of the two RC layers, $f_{\text{in}}$ in Eq. (4). MZM 1 is
driven by an arbitrary waveform generator (AWG 1) which supplies the input
signal $u_{n}^{(1)}$; MZM 2 can be driven both by a second arbitrary waveform
generator (AWG 2) or by the output of a photodiode (PD 2), as described below.
The two modulated signals are merged together in a 50/50 fiber coupler (C1)
and then injected into an Erbium-doped-fiber amplifier (EDFA 1). EDFA 1 raises
the total power to $9$ dBm, equally distributed between the two signals. After
the amplification, the two signals are injected in a phase modulator (PM 1).
PM 1 is driven by a sinusoidal radio-frequency signal (frequency
$\Omega\approx 17$ GHz, power P1$\approx 30$ dBm). The radio-frequency signal
is generated by an RF clock (RF source) and amplified by an RF amplifier (RF
AMP 1). The phase modulation provided by PM 1 generates two frequency combs
centered in $\lambda_{1}$ and $\lambda_{2}$ (Fig. 3). The spacing of the comb
lines is equal to $\Omega$ and the number of lines depends on P1. In our
implementation, PM 1 provides approximately $20$ usable comb lines per comb,
i.e., $20$ neurons. The two combs constitute the input stimuli for the two
reservoir networks. The amplitude of each line determines how strongly the
input signal is coupled to the particular neuron encoded in that line. Hence,
the distribution of (complex) amplitudes among the comb lines defines the two
vectors of input-to-reservoir weights, $\mathbf{W}_{\text{in}}^{\text{(1)}}$
and $\mathbf{W}_{\text{in}}^{\text{(2)}}$. The two frequency combs are
injected in a fiber loop through a 30/70 coupler (C2). The fiber loop is 15
meters long, corresponding to a roundtrip frequency of approximately $20$ MHz.
All the input signals are synchronized with the roundtrip time of the loop, in
such a way that each timestep of the input signals entirely fills the loop.
Hence, the processing frequency of our system is fixed by the cavity length
and is approximately $20$ MHz. The fiber loop contains a second phase
modulator (PM 2) and an optical amplifier (EDFA 2). PM 2 is driven by a signal
generated by the same RF source but undergoes a different amplification (RF
AMP 2), hence it has the same frequency but a different power P2$\approx$20
dBm as the RF signal supplied to PM 1. The phase modulation provided by PM 2
creates frequency interference among the lines of the same comb, thus
implementing the (complex-weighted) connectivity among the neurons of the same
reservoir. EDFA 2 compensates for the losses in the loop. The transformation
of the combs over a roundtrip, including the effects of phase modulation,
amplification and dispersion (which acts differently on each comb line/neuron)
define the matrices $\mathbf{W}_{\text{}}^{\text{(1)}}$ and
$\mathbf{W}_{\text{}}^{\text{(2)}}$. The amplitudes of the two combs at each
roundtrip ${n}$ provide the states of the two reservoirs
$\mathbf{x}_{n}^{\text{(1)}}$ and $\mathbf{x}_{n}^{\text{(2)}}$.

Part of the circulating radiation is extracted by a 20/80 fiber coupler (C3),
amplified by EDFA 3 and directed to the readout circuit. The readout consists
of a multi-channel programmable spectral filter (PSF, Coherent II-VI
Waveshaper) and two photodiodes (PD 1 and PD 2), measuring each of the two PSF
outputs. The first PSF channel, connected to PD 1, is employed to measure the
evolution of both reservoirs. The measurement procedure consists of selecting
a single comb line using a band-pass filter and recording via PD 1 the
intensity of that comb line. At the end of the procedure, the intensities of
all the comb lines, i.e. the norm square of the components of vectors
$\mathbf{x}_{n}^{\text{(1)}}$ and $\mathbf{x}_{n}^{\text{(2)}}$, are recorded
on a computer. The output of the reservoir is then obtained by multiplying
these intensities by the output weights
$|\mathbf{W}_{\text{out}}^{\text{(\text{T+})}}|^{2}$ and
$-|\mathbf{W}_{\text{out}}^{\text{(\text{T-})}}|^{2}$. Note that this
operation could be realised by using two channels of the PSF to which one
applies weights $\mathbf{W}_{\text{out}}^{\text{(\text{T+})}}$ and
$\mathbf{W}_{\text{out}}^{\text{(\text{T-})}}$ respectively, sending the two
outputs to two photodiodes, and taking the difference of the resulting
currents, as described in [6].

##### Operation modes.

We use two operation modes: “deep” and “independent”.

In deep-RC mode, the second channel of the programmable spectral filter is
configured to select and transmit only the comb centered on $\lambda_{1}$
after having applied an attenuation mask
$\mathbf{W}_{\text{out}}^{\text{(1)}}$. Consequently, PD 2 measures the signal
$u_{n}^{(2)}=\left|\mathbf{W}_{\text{out}}^{\text{(1)}}\cdot\mathbf{x}_{n}^{\text{(1)}}\right|^{2}$.
The output of PD 2 drives MZM 2, and thus constitutes the input of the second
RC. In this configuration the system is a two layer deep-RC, as described in
subsection 2.1.2.

In independent mode, the two RC computations are decoupled by driving MZM 2
through a second, independent, arbitrary waveform generator AWG 2 (the second
channel of PSF and PD 2 are deactivated). The two RC computations do not
interact with each other and are carried out independently.

The selection of the computation mode, deep or independent, is made by
flipping an electric switch that selects whether MZM 2 is driven by PD 2 or by
AWG 2, as illustrated in Fig. 2.

##### Stabilization.

The experimental setup is sensitive to acoustical noise and thermal drift. To
limit these effects, the optical loop, including PM 2 and EDFA 2, is mounted
inside an insulated box on an optical table. Furthermore, two PID controllers
piezo-tune the emission wavelengths of the two laser sources in order to fix
the operating condition to a certain point in the loop transfer function. The
two PID controllers are fed by the intensity of the reflection of each of the
two combs at the entrance of the loop, on the coupler C2. This requires two
auxiliary photodiodes and spectral filters (not represented in Fig. 2).

###  2.3 Benchmark tasks

We selected two benchmark tasks, the first consisting of the prediction of the
evolution of a chaotic time series, and the second one consisting of the
compensation of the distortion produced in a nonlinear communication channel.

The time series prediction task is based on the infrared laser dataset of the
Santa Fe Time Series Competition [36]. The time series $u_{t}$ is supplied as
input, and the task consists of producing $u_{t-\tau}$, with
$-5\leq\tau\leq+5$. Note that when the timeshift $\tau$ is negative, the task
consists of remembering the past, while when the timeshift $\tau$ is positive,
the task consists of predicting the future. The accuracy is expressed in terms
of Normalized Mean Square Error (NMSE) between the target signal and the
produced output. When running this benchmark, the training set is composed of
6000 timesteps, and the testing set is composed of 2500 timesteps (this is a
standard 70%-30% repartition). We discard the first 500 timesteps of reservoir
output to avoid operating in a transient regime.

The nonlinear compensation task was first used in the RC community in [14]. A
random signal composed of four different symbols is propagated along a
simulated channel exhibiting nonlinearity, noise, and memory about past
inputs. The task consists in reconstructing the original input given the
channel output. Performance is evaluated for different Signal-to-Noise Ratios
(SNR) in the range $[8\text{dB},\text{ }32\text{dB}]$. The results are
expressed in terms of Symbol Error Rate (SER), i.e., the ratio of wrongly
reconstructed output symbols over the total number of transmitted symbols.
When running this benchmark, the training set is composed of 14000 timesteps,
and the testing set is composed of 30000 timesteps. We discard the first 1000
timesteps of reservoir output to avoid operating in a transient regime. Note
that, differently from the chaotic time series benchmark which is based on a
limited dataset, in this case the training dataset can be easily generated at
the moment. This is why we employed a larger amount of datapoints for the
initial wash-out and the testing.

Every benchmark result has been validated through a 100-steps cross-
validation, meaning that the points belonging to the train and test datasets
have been selected at random for 100 times and results have been averaged.

###  2.4 Tested configurations

(a)

(b)

(c)

Figure 4: The three tested configurations for the two independent RCs. “
Reservoir $\lambda_{1}$” is encoded in the frequency comb centered around
$\lambda_{1}$, while “reservoir $\lambda_{2}$” is encoded in the frequency
comb centered around $\lambda_{2}$. Both reservoirs are executed on the same
photonic substrate. (a) Shallow-RC: one of the two reservoirs performs the
benchmark task as a traditional RC, while the other reservoir processes a
different time series in parallel. (b) Parallel-RC: both reservoirs process
the same input time series, but their dynamics are decoupled from each other.
A single output layer is trained, which combines signals from both reservoirs.
(c) Deep-RC: the two reservoirs constitute the two layers of a deep-RC.

Our photonic system supports two RCs that operate simultaneously, either
independently from each other, or connected in series. We evaluated the
performance of three different configurations, described in Fig. 1. First, a
“shallow-RC” configuration (Fig. 4(a)) in which only one of the two
independent RCs executes the benchmark task, constituting a “traditional” RC
as described in Sec. 2.1.1. In this configuration, the second RC processes a
different, not evaluated, computation, with the purpose of simulating a
parallel computation scenario where two different tasks are performed at the
same time. Second, a “parallel-RC” configuration (Fig. 4(b)) where both
independent RCs execute the same task in an uncorrelated way, and a single
output layer is connected to both reservoirs. This constitutes a “non-deep”
way of using the full computational capabilities of the system on a single
task. Third, a “deep-RC” configuration (Fig. 4(c)) where the two independent
RCs are connected in series as described in Sec. 2.1.2.

In addition, in the deep-RC configuration, we used two different methods to
tune the weights $\mathbf{W}_{\text{out}}^{\text{(1)}}$, i.e. the attenuations
applied by the PSF, that determine the connection from the first RC layer to
the second one. In the first, simplest, approach, we apply the same
attenuation to all comb lines, corresponding to
$\mathbf{W}_{\text{out}}^{\text{(1)}}=\textrm{diag}(\alpha)$, and we optimize
the overall attenuation $\alpha^{2}$ by sweeping it in the range
$[-20\text{dB},\text{ }0\text{dB}]$. In the second approach, we optimise all
the coefficients of $\mathbf{W}_{\text{out}}^{\text{(1)}}$ by using the
Covariance Matrix Adaptation Evolution Strategy (CMA-ES) optimisation
algorithm [12].

Finally, in order to improve the reservoir computing performance, we tune the
comb line spacing $\Omega$ to the best-performing value for each task. The
fiber loop constitutes a spectral interferometer and exhibits, due to
dispersion in the fiber, a complex behavior strongly dependent on $\Omega$.
This is illustrated in Fig. 5, where the performance in shallow-RC
configuration on the two tasks is plotted as a function of $\Omega$ for both
reservoirs.

Figure 5: Performance of the reservoir computers in shallow-RC configuration
as a function of $\Omega$ on the channel equalization task (top) and the
Santa-Fe time series prediction task for prediction 1 timestep ahead (bottom)
. The complex dependence on $\Omega$ is due to the dispersion in the optical
fiber. The dispersion is also the reason why the dependence on $\Omega$ is
different for RC-1 and RC-2, as they use frequency combs centered on different
wavelengths. (As these plots are time-consuming to obtain, a reduced number of
comb lines $N=14$ was used).

##  3 Results

(a)

(b)

Figure 6: Experimental results for the three operation modes (shallow-RC,
parallel-RC, and deep-RC) on the two selected benchmark tasks: nonlinear
channel equalization (a) and chaotic time series prediction (b). Deep-RC
results are shown for both optimization methods presented in the text (uniform
optimized attenuation $\alpha$ and CMA-ES). Error-bars represent the score
standard deviation measured in cross-validation phase. Results in (a) are
expressed as symbol error rate (SER) vs. signal to noise ratio (SNR). Results
in (b) are expressed as normalized mean square error (NMSE) vs. shift of the
target time series with respect to the input one. When the shift is positive,
the task consists in predicting the future; when the shift is zero, the task
consists in reproducing the present input; when the shift is negative, the
task consists in reproducing the past.

Results of the benchmark tests are reported in Fig. 6 for the three operation
modes: shallow-RC, parallel-RC, and deep-RC. In this figure, the deep-RC
results are shown, for both the optimization techniques described in Sec. 2.4.
Nonlinear channel equalization results (Fig. 6(a)) show the expected decrease
of Symbol Error Rate (SER) when the Signal-to-Noise Ratio (SNR) increases,
since the addition of noise increases the difficulty of the task, eventually
making it impossible to revert the channel distortion. For high-SNR values,
both shallow-RC and parallel-RC SER scores saturate, while deep-RC SER score
maintains an exponential decay for increasing SNR values. For every SNR value,
deep-RC always performed better, followed by the parallel-RC and finally by
shallow-RC. A similar behavior is found in the results of the chaotic time
series prediction task (Fig. 6(b)).

Two trends are clearly visible from Fig. 6. First, parallel-RC systematically
outperforms shallow-RC. Indeed since the two parallel RCs perform different
computations (as is evident from Fig. 5), using both reservoirs in parallel
should perform at least as well as using a single reservoir. Second, deep-RC
outperforms parallel-RC in every test we conducted. Both configurations
exploit the same number of neurons and are differentiated only by the
topology, thus we conclude that the serial configuration in deep-RC really
boosts RC performance.

We observe that the two optimization techniques for the inter-layer connection
perform comparably, with the simpler algorithm sometimes outperforming the
CMA-ES algorithm. We identified two reasons for this behavior. First, the CMA-
ES algorithm could get stuck in local minima. Second, the search of the
optimal set of weights could be affected negatively by slow drifts in the
operating conditions of the deep-RC.

##  4 Conclusion

We presented a fully analog photonic implementation of a deep reservoir
computer. The connection between the two layers is performed in the analog
domain with no processing or storing on a digital computer. The presented
implementation also allows for two independent RC computations to be executed
at the same time. We found that the deep-RC configuration, obtained by
connecting in series the two RCs, performs better than a parallel-RC
configuration, where the two RCs are connected to the same output, and process
the same input data, but do not interact during computation.

The reported experiment has only two layers, but deeper schemes are in
principle possible. New layers can be added to the deep-RC by using more than
two lasers, provided that the generated combs do not overlap each other. The C
band could host $10$ parallel computations (considering combs $3$ nm wide, see
Fig. 3). These $10$ parallel computations could be employed to constitute a
single $10$-layers deep-RC, or even multiple deep-RC running in parallel, each
one composed of fewer layers. On the other hand, broader combs would be able
to encode more neurons in each reservoir. Thus, a balance between the number
of layers and the number of neurons per layer has to be searched for. In any
case, integrating (partially or entirely) the experiment, as proposed in [17],
could be a route to scaling up the system while simplifying its stabilization.

Although we already explored two strategies for optimizing the interconnection
among the two deep-RC layers, many ideas are still to be tested (see e.g.,
[10, 38, 25]) and could be the object of further investigation.

In summary, developing deep architectures for neuromorphic photonic computing
is a highly promising avenue for increasing both the complexity of the tasks
that can be solved and the system performance. However, the presence of
analog-to-digital or digital-to-analog converters strongly affects power
consumption and footprint, hence it is to be avoided. We have demonstrated
that this is possible for photonic deep reservoir computing.

## Funding

The authors acknowledge financial support from the European Union’s Horizon
2020 research and innovation program under the Marie Skłodowska-Curie grant
agreement 860830 (POST DIGITAL), and from the FWO and F.R.S.-FNRS Excellence
of Science (EOS) programme grant 40007536.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available
at this time but may be obtained from the authors upon reasonable request.

## References

  * [1] Lennert Appeltant et al.  “Information processing using a single dynamical node as complex system”  In _Nature communications_ 2.1 Nature Publishing Group, 2011, pp. 1–6
  * [2] Farshid Ashtiani, Alexander J Geers and Firooz Aflatouni  “An on-chip photonic deep neural network for image classification”  In _Nature_ 606.7914 Nature Publishing Group UK London, 2022, pp. 501–506
  * [3] Minkyung Baek et al.  “Accurate prediction of protein structures and interactions using a three-track neural network”  In _Science_ 373.6557, 2021, pp. 871–876  DOI: 10.1126/science.abj8754
  * [4] Y. Bengio  “Learning Deep Architectures for AI”  In _Foundations and Trends in Machine Learning_ 2.1 Now Publishers, 2009, pp. 1–127  DOI: 10.1561/2200000006
  * [5] “Bing Chat” Accessed: March 2023, https://www.bing.com
  * [6] Lorenz Butschek et al.  “Photonic reservoir computer based on frequency multiplexing”  In _Optics Letters_ 47.4 Optica Publishing Group, 2022, pp. 782  DOI: 10.1364/ol.451087
  * [7] Johannes Feldmann et al.  “All-optical spiking neurosynaptic networks with self-learning capabilities”  In _Nature_ 569.7755 Nature Publishing Group UK London, 2019, pp. 208–214
  * [8] Johannes Feldmann et al.  “Parallel convolutional processing using an integrated photonic tensor core”  In _Nature_ 589.7840 Nature Publishing Group UK London, 2021, pp. 52–58
  * [9] Matthias Freiberger et al.  “Improving time series recognition and prediction with networks and ensembles of passive photonic reservoirs”  In _IEEE Journal of Selected Topics in Quantum Electronics_ 26.1 IEEE, 2019, pp. 1–11
  * [10] Claudio Gallicchio, Alessio Micheli and Luca Pedrelli  “Deep reservoir computing: A critical experimental analysis”  In _Neurocomputing_ 268 Elsevier BV, 2017, pp. 87–99  DOI: 10.1016/j.neucom.2016.12.089
  * [11] Ryan Hamerly et al.  “Large-scale optical neural networks based on photoelectric multiplication”  In _Physical Review X_ 9.2 APS, 2019, pp. 021032
  * [12] Nikolaus Hansen  “The CMA evolution strategy: a comparing review”  In _Towards a new evolutionary computation: Advances in the estimation of distribution algorithms_ Springer, 2006, pp. 75–102
  * [13] Guang-Bin Huang, Qin-Yu Zhu and Chee-Kheong Siew  “Extreme learning machine: theory and applications”  In _Neurocomputing_ 70.1-3 Elsevier, 2006, pp. 489–501
  * [14] Herbert Jaeger and Harald Haas  “Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication”  In _Science_ 304.5667 American Association for the Advancement of Science (AAAS), 2004, pp. 78–80  DOI: 10.1126/science.1091277
  * [15] Aashu Jha et al.  “Photonic spiking neural networks and graphene-on-silicon spiking neurons”  In _Journal of Lightwave Technology_ 40.9 IEEE, 2022, pp. 2901–2914
  * [16] John Jumper et al.  “Highly accurate protein structure prediction with AlphaFold”  In _Nature_ 596.7873 Nature Publishing Group UK London, 2021, pp. 583–589
  * [17] Wosen Kassa et al.  “Towards integrated parallel photonic reservoir computing based on frequency multiplexing”  In _Neuro-inspired Photonic Computing_ SPIE, 2018  DOI: 10.1117/12.2306176
  * [18] Laurent Larger et al.  “High-speed photonic reservoir computing using a time-delay-based architecture: Million words per second classification”  In _Physical Review X_ 7.1 APS, 2017, pp. 011015
  * [19] Yann LeCun, Yoshua Bengio and Geoffrey Hinton  “Deep learning”  In _Nature_ 521.7553 Springer ScienceBusiness Media LLC, 2015, pp. 436–444  DOI: 10.1038/nature14539
  * [20] Bao-De Lin et al.  “Deep time-delay reservoir computing with cascading injection-locked lasers”  In _IEEE Journal of Selected Topics in Quantum Electronics_ 29.6: Photonic Signal Processing IEEE, 2022, pp. 1–8
  * [21] Xing Lin et al.  “All-optical machine learning using diffractive deep neural networks”  In _Science_ 361.6406 American Association for the Advancement of Science, 2018, pp. 1004–1008
  * [22] Antoine Liutkus et al.  “Imaging With Nature: Compressive Imaging Using a Multiply Scattering Medium”  In _Scientific Reports_ 4 Nature Publishing Group UK London, 2014, pp. 5552
  * [23] Danijela Marković, Alice Mizrahi, Damien Querlioz and Julie Grollier  “Physics for neuromorphic computing”  In _Nature Reviews Physics_ 2.9 Nature Publishing Group UK London, 2020, pp. 499–510
  * [24] Mitsumasa Nakajima, Kenji Tanaka and Toshikazu Hashimoto  “Scalable reservoir computing on coherent linear photonic processor”  In _Communications Physics_ 4.1 Nature Publishing Group UK London, 2021, pp. 20
  * [25] Mitsumasa Nakajima et al.  “Physical deep learning with biologically inspired training method: gradient-free approach for physical hardware”  In _Nature Communications_ 13.1 Springer ScienceBusiness Media LLC, 2022  DOI: 10.1038/s41467-022-35216-2
  * [26] “OpenAI” Accessed: March 2023, https://openai.com/about
  * [27] Mushegh Rafayelyan et al.  “Large-scale optical reservoir computing for spatiotemporal chaotic systems prediction”  In _Physical Review X_ 10.4 APS, 2020, pp. 041037
  * [28] Sandeep Ravidran  “Five ways deep learning has transformed image analysis”  In _Nature_ 609, 2022, pp. 864–866
  * [29] A. Saade et al.  “Random projections through multiple optical scattering: Approximating Kernels at the speed of light”  In _2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2016, pp. 6215–6219  DOI: 10.1109/ICASSP.2016.7472872
  * [30] Yichen Shen et al.  “Deep learning with coherent nanophotonic circuits”  In _Nature photonics_ 11.7 Nature Publishing Group UK London, 2017, pp. 441–446
  * [31] David Silver et al.  “A general reinforcement learning algorithm that masters chess, shogi, and Go through self-play”  In _Science_ 362.6419 American Association for the Advancement of Science, 2018, pp. 1140–1144
  * [32] Satoshi Sunada, Kazutaka Kanno and Atsushi Uchida  “Using multidimensional speckle dynamics for high-speed, large-scale, parallel photonic computing”  In _Optics Express_ 28.21 Optical Society of America, 2020, pp. 30349–30361
  * [33] Gouhei Tanaka et al.  “Recent advances in physical reservoir computing: A review”  In _Neural Networks_ 115 Elsevier BV, 2019, pp. 100–123  DOI: 10.1016/j.neunet.2019.03.005
  * [34] Fabian Triefenbach, Azarakhsh Jalalvand, Benjamin Schrauwen and Jean-Pierre Martens  “Phoneme recognition with large hierarchical reservoirs”  In _Advances in neural information processing systems_ 23, 2010
  * [35] Quentin Vinckier et al.  “High-performance photonic reservoir computer based on a coherently driven passive cavity”  In _Optica_ 2.5 Optical Society of America, 2015, pp. 438–446
  * [36] Carl-Otto Weiss, Udo Hübner, Neal Broadus Abraham and Dingyuan Tang  “Lorenz-like chaos in NH3-FIR lasers”  In _Infrared Physics & Technology_ 36.1 Elsevier BV, 1995, pp. 489–512  DOI: 10.1016/1350-4495(94)00088-3
  * [37] Gordon Wetzstein et al.  “Inference in artificial intelligence with deep optics and photonics”  In _Nature_ 588.7836 Nature Publishing Group UK London, 2020, pp. 39–47
  * [38] Logan G Wright et al.  “Deep physical neural networks trained with backpropagation”  In _Nature_ 601.7894 Nature Publishing Group UK London, 2022, pp. 549–555
  * [39] Xingyuan Xu et al.  “11 TOPS photonic convolutional accelerator for optical neural networks”  In _Nature_ 589.7840 Nature Publishing Group UK London, 2021, pp. 44–51
  * [40] Tiankuang Zhou et al.  “Large-scale neuromorphic optoelectronic computing with a reconfigurable diffractive processing unit”  In _Nature Photonics_ 15.5 Nature Publishing Group UK London, 2021, pp. 367–373