# Weak Kerr Nonlinearity Boosts the Performance of Frequency-Multiplexed
Photonic Extreme Learning Machines: A Multifaceted Approach

Marina Zajnulina  1 Alessandro Lupo  2,* and Serge Massar2  1 Multitel
Innovation Centre, Rue Pierre et Marie Curie 2, 7000 Mons, Belgium
2Laboratoire d’Information Quantique CP224, Université libre de Bruxelles
(ULB), Av. F. D. Roosevelt 50, 1050 Bruxelles, Belgium
<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

We provide a theoretical, numerical, and experimental investigation of the
Kerr nonlinearity impact on the performance of a frequency-multiplexed Extreme
Learning Machine (ELM). In such ELM, the neuron signals are encoded in the
lines of a frequency comb. The Kerr nonlinearity facilitates the randomized
neuron connections allowing for efficient information mixing. A programmable
spectral filter applies the output weights. The system operates in a
continuous-wave regime. Even at low input peak powers, the resulting weak Kerr
nonlinearity is sufficient to significantly boost the performance on several
tasks. This boost already arises when one uses only the very small Kerr
nonlinearity present in a 20-meter long erbium-doped fiber amplifier. In
contrast, a subsequent propagation in 540 meters of a single-mode fiber
improves the performance only slightly, whereas additional information mixing
with a phase modulator does not result in a further improvement at all. We
introduce a model to show that, in frequency-multiplexed ELMs, the Kerr
nonlinearity mixes information via four-wave mixing, rather than via self- or
cross-phase modulation. At low powers, this effect is quartic in the comb-line
amplitudes. Numerical simulations validate our experimental results and
interpretation.

##  1 Introduction

Artificial neural networks (ANNs) are brain-inspired algorithms that process
information by propagating it through a series of weighted connections and
nonlinear activation functions. ANNs are usually implemented in electronic
digital computers, mainly through matrix multiplications with a subsequent
application of nonlinear functions. Due to the separation between computing
and memory units, the classical von Neumann computer architecture suffers a
severe performance bottleneck when executing algorithms that massively depend
on memory access, such as ANNs [1]. For this reason, ANN tasks are usually
executed on auxiliary hardware (GPUs or AI accelerators)  [2]. In these
devices, the information flows through an artificial structure of hard-wired
networks and/or logic gates. In opposition, the biological brain exploits the
complex dynamics of its organic material. This inspires alternative approaches
to computing that harvest the dynamics of the physical substrate they are
built on, leading to the implementation of analog “physical computers” [3, 4].

Extreme Learning Machines (ELMs) are feed-forward ANNs with fixed, i.e. not
trained, connections between the input and hidden layers, while only the
output-layer connections undergo a training procedure to find optimal
weights[5, 6]. The output-layer training can be done using non-iterative
algorithms such as linear regression, which is by far less energy- and time-
consuming than the backpropagation algorithm usually employed for training of
fully connected networks. It has been proven that training only the output
connections is enough to guarantee the universal approximation property of
ELMs [7], although this is only a necessary (but not sufficient) condition for
good performance. Since there is no requirement for internal-connections
tuning, ELMs are particularly suitable for implementation on physical
substrates. In such physical implementations, the substrate nonlinearly mixes
the input signals and therefore plays the role of connectivity between the
input and hidden layers. Thus, Ref. [8] reports a theoretical study of the
computing capability of nonlinear waves such as rogue waves, dispersive
shocks, and solitons.

Photonic systems constitute an interesting substrate for analog computation.
They are not only able to provide different types of nonlinearity for
information processing, but also exhibit multiple degrees of freedom for
information encoding, as well as the potential for low-loss, ultra-fast, and
high-bandwidth data processing and throughput[9]. For example, photonic ELMs
have been implemented based on propagation through scattering media [10, 11]
and multimode fibers [12], as well as in delay-based schemes [13], and, in
general, in any system whose dynamics is described by the Nonlinear
Schrödinger equation (NLS) [8].

Using the frequency of light as a degree of freedom is a natural solution for
processing large numbers of signals in parallel. Frequency-multiplexed modes
can be generated and manipulated by employing reliable and efficient optical
components, such as phase and intensity modulators, programmable spectral
filters, and other devices for wavelength-division multiplexing (WDM). All
these components can be integrated on a photonic chip and are already
massively employed for telecommunication. Recently, various frequency-
multiplexing-based schemes have been used to implement physical ANNs [14, 15,
16, 17, 18, 19, 20, 21, 22, 23], where, depending on the experiment, the
frequency-multiplexed information is manipulated by wavelength-dependent
couplers, or exploiting the dispersion or optical nonlinearities of the
material.

In our previous work [17, 18], we introduced a photonic ELM that uses
frequency multiplexing to encode the input information in the lines of an
optical frequency comb. The information mixing in the input-to-hidden layer
connections was realized via frequency interference induced by periodical
phase modulation, giving rise to a transformation linear in the electrical
field amplitudes, i.e. we had a linear mixing of input signals. The required
nonlinearity was quadratic. It was provided by a photodiode that measured the
intensities of the comb lines.

Here, we exploit the Kerr nonlinearity as a fundamentally different
information processing mechanism. It arises in single-mode optical fibers
(erbium-doped and standard telecom fibers) deployed in the experiment. Thus,
we use a low-power continuous wave (CW) operation in an approximately
$20$-meter-long erbium-doped fiber optionally followed by a $540$-meter-long
standard telecom fiber.

We demonstrate our system on four benchmark tasks, consisting of
classification and regression problems. We also compare our system to
software-based ELMs with the same number of nodes, obtaining similar
prediction accuracy. Finally, we present an analytical model that shows that
the information mixing in our system occurs due to four-wave mixing (FWM);
dispersion effects as well as self- and cross-phase modulation play no role.
In the low-power regime, the Kerr effect introduces a quartic nonlinearity (as
opposed to quadratic nonlinearity from our previous work).

We find that using the Kerr nonlinearity generated in a piece of only
20-meter-long erbium-doped fiber amplifier (EDFA) is sufficient to
significantly boost the performance of our ELM on several tasks. Subsequent
propagation in 540 meters of a single-mode fiber improves the performance only
slightly. Linear information mixing with a following phase modulator does not
result in a further improvement at all. As the measured EDFA output is only in
the rage of 0.5W, we deal with a weak overall Kerr nonlinearity in the system.
This means that a weak Kerr nonlinearity is sufficient for effective and
efficient information processing.

The sufficiency of a weak Kerr nonlinearity is a counterintuitive result
opposing the general idea that a high nonlinearity is required for photonic
neuromorphic computing (cf. [8]). In fact, Refs. [23] and [21] report using
highly nonlinear fibres (HNLF) and femto-second (pulsed) lasers as sources to
facilitate information processing. Highly nonlinear regimes require high peak
power pulsed lasers and/or special fibers, effectively reducing the
applicability of the approach. Using a weak Kerr nonlinearity makes the system
easier to implement. We note that a successful usage of a weak Kerr
nonlinearity was also studied in Ref. [24] in the context of Reservoir
Computing and in Refs. [25, 26] in the context of Hopfield networks, showing
the generalization potential of this concept. In addition, fiber-based schemes
that deploy weak Kerr nonlinearities would be easier to transfer on
semiconductor chips than the schemes that operate in highly nonlinear regimes,
as they require less optical power and yield less thermal dissipation.

The paper is structured as follows. Sec. 2 describes the principles of a
generic ELM algorithm, the impact of the Kerr nonlinearity and frequency-
domain interference; Sec. 3 describes the experimental setup, as well as its
operating procedures; Sec. 4 contains the experimental results and their
discussion; Sec. 5 contains the conclusions. App. A depicts the low-power
analytical model for information mixing driven by optical Kerr effect. App. B
reports a summary of the benchmarking results in every tested configuration.

##  2 Principles

###  2.1 Extreme Learning Machine Algorithm

An Extreme Learning Machine (ELM) is a feed-forward neural network where only
the output weights are trained, while the other connections are fixed and
remain untrained [5, 6]. We define $N^{\mathit{in}},$ $N^{\mathit{hidden}},$
and $N^{\mathit{out}}$ as the numbers of the input, hidden, and output
neurons, respectively. $\mathbf{u}$, $\mathbf{h}$, and $\mathbf{y}$ are the
column vectors containing the corresponding neuron values. In general, a
single-hidden–layer ELM is described by the following two equations:

$\displaystyle\mathbf{h}$
$\displaystyle=f_{\mathit{hidden}}\left(\mathbf{W}^{\mathit{in}}\cdot
f_{\mathit{in}}\left(\mathbf{u}\right)\right),$ (1) $\displaystyle\mathbf{y}$
$\displaystyle=\mathbf{W}^{\mathit{out}}\cdot\mathbf{h},$ (2)

where $f_{\mathit{in}}$ and $f_{\mathit{hidden}}$ are the nonlinear activation
functions of the input and hidden layers (they act elementwise),
$\mathbf{W}^{\mathit{in}}$ is a $N_{\mathit{hidden}}\times N_{\mathit{in}}$
matrix of the input-to-hidden layer weights, and $\mathbf{W}^{\mathit{out}}$
is a $N^{\mathit{output}}\times N^{\mathit{hidden}}$ matrix of the hidden-to-
output layer weights. The matrix $\mathbf{W}^{\mathit{in}}$ is not trained,
but rather selected at random. $\mathbf{W}^{\mathit{out}}$ is trained
employing ridge regression to achieve an output $\mathbf{y}$ that approximates
the true value.

In physical implementations of ELMs, the transformation
$\mathbf{h}(\mathbf{u})$ between the input and the hidden layer is fixed by
the physical system and may therefore take a form different from Eq. (1).
However it is important for good operation that $\mathbf{h}(\mathbf{u})$ is a
nonlinear mapping from the $N^{\mathit{in}}$ input nodes to a larger number of
hidden nodes $N^{\mathit{hidden}}\gg N^{\mathit{in}}$. In the context of ANNs,
this mapping corresponds to the mapping of the input to a higher-dimensional
space which generally allows for a better data separability.

###  2.2 Extreme Learning Machine based on Frequency Multiplexing

#### General Framework.

As in our previous publications (cf. [17, 18]), we use an optical frequency
comb with central frequency $\omega_{0}$ and line spacing $\Omega$ to encode
and process the input data. The k-th line in the comb reads accordingly as
$\omega_{k}=(\omega_{0}+k\Omega)$, where $k$ is a positive or negative integer
representing frequencies higher or lower than $\omega_{0}$.

We encode the input data $\mathbf{u}$ (cf. Eq. 1) in the comb-line amplitudes
$x^{\textit{in}}_{k}$ by applying attenuations or phase delays with
corresponding values. More details on data encoding can be found in Sec. 3.2.
With encoded input data, the time-domain representation of the comb reads as:

$E^{\textit{in}}(t)=\sum_{k}x^{\textit{in}}_{k}(\mathbf{u})e^{-i\left(\omega_{0}+k\Omega\right)t}\
,$ (3)

where $t$ represents the time, $k$ is the comb line index with $k=0$
corresponding to the central line, and $x^{\textit{in}}_{k}$ is the k-th line
complex amplitude. Since the frequency comb has a limited spectral range, the
elements of the sum in Eq. (3) vanish for a sufficiently large $|k|$.

After propagation through the system, the electric field is transformed into
the output field:

$E^{\textit{out}}(t)=\sum_{k}x^{\textit{out}}_{k}(\mathbf{u})e^{-i\left(\omega_{0}+k\Omega\right)t}\
,$ (4)

where the output amplitudes are functions of the input amplitudes:

${x}^{\textit{out}}_{k}=f_{k}(\mathbf{x}^{\textit{in}})\ .$ (5)

In our previous work, the transformation Eq. (5) was linear and realized by a
phase modulator [17]. Here, the transformation is facilitated by Kerr
nonlinearity as it will be shown in the further course of our study.

We obtain the ELM output by taking a linear combination of the comb line
intensities:

$y_{\textit{j}}=\sum_{k}W^{\textit{out}}_{jk}|x^{\textit{out}}_{k}|^{2}\ .$
(6)

where $\mathbf{W}^{\textit{out}}$ is trained by ridge regression. Note that
the operation Eq. (6) can be realised optically, as reported in [17] and
described below. In summary, the vector of output comb intensities has the
meaning of a hidden layer $\mathbf{h}=|\mathbf{x^{\textit{out}}}|^{2}$, with
Eqs. (3, 4, 5, 6) corresponding to Eqs. (1, 2).

#### Fiber Kerr Nonlinearity for Information Mixing.

Figure 1:  Evolution of five most powerful comb-line intensities normalized by
comb total power with fiber parameters taken from Sec. 3.1.
(a) Normalized comb-line intensities as a function of fiber propagation length
for a fixed input power of $0.5$ W. Solid lines are obtained by integration of
full Nonlinear Schrödinger equation via split-step Fourier method (SSFM).
Dashed lines show intensities obtained via integration of comb-amplitude
equation derived explicitly for our Extreme Learning Machine (Eq. (23) in App.
A). The two methods should be equivalent and the slight diffference is
probably due to numerical precision.
(b) Same as (a), but dashed lines are produced by integration of Eq. (31)
obtained via perturbing the Nonlinear Schrödinger equation and taking the
terms in the order of $\gamma$ (App. A). The figure shows that the linearized
approximation of Eq. (31) is valid for approximately $100$ meters of
propagation.
(c) Normalized comb-line intensities as a function of a function of the input
power for a fixed fiber propagation length of $L=540$ m obtained via
integration of the Nonlinear Schrödinger equation.

The input transforms to the output as the comb propagates through an optical
fiber of length $L$ (cf. Eq. (5)). During this propagation, the slowly varying
electric-field amplitude $E(z,t)$ formulated in the co-moving frame undergoes
linear (dispersive) and nonlinear (Kerr-effect-based) transformations
described by NLS [8, 27]:

$\frac{\partial E}{\partial
z}=-i\frac{\beta_{2}}{2}\frac{\partial^{2}E}{\partial
t^{2}}+i\gamma|E|^{2}E-\frac{\alpha}{2}{E},$ (7)

where $E(z,t)$ is the electric-field amplitude at position $z$ and time $t$,
$\beta_{2}$ is the group-velocity dispersion, $\gamma$ is the nonlinear
coefficient that describes the intensity-dependent Kerr effect, and $\alpha$
are the optical losses (with $\alpha<0$ corresponding to a gain, as in the
case of an EDFA).

The electric field at the output of the optical fiber reads as:

$E(L,t)=\sum_{k}x^{\textit{Kerr}}_{\textit{k}}e^{-i(\omega_{0}+k\Omega)t}.$
(8)

where $x^{\textit{Kerr}}_{k}$ are comb-line amplitudes transformed by Kerr
nonlinearity. The corresponding $j$-th output of the ELM takes the form:

$y_{\textit{j}}=\sum_{\textit{k}}W^{out}_{\textit{jk}}|x^{\textit{Kerr}}_{\textit{k}}|^{2}\
.$ (9)

Eq. 7 is not integrable. Therefore, we cannot provide an analytical expression
for the dependence between the input $(x^{\textit{in}}_{k})$ and output
$(x^{\textit{Kerr}}_{\textit{k}})$ comb lines. However, we derived a model for
the Kerr-based comb-line transformation assuming a weak Kerr nonlinearity,
i.e. a nonlinearity of the first order in $\gamma$. More details on this model
can be found in App. A.

Our model shows that

  * •

Kerr-nonlinearity-driven FWM is the only process that is responsible for
information mixing in our ELM. It affects the intensities of the comb lines,
which is detected by ELM’s readout layer. In contrast, self- and cross-phase
modulation do not contribute to information mixing as these effects act on the
phases of $x^{\textit{Kerr}}_{\textit{k}}$ and the phase information is lost
when the comb line intensities are read out.

  * •

The change in the comb line intensities, i.e. the information mixing, depends
quartically on comb line amplitudes If we further suppose that the porpagation
distance $L$ is very small, so that dispersion and loss can be neglected, then
the output intensities take the form

$|x^{\textit{Kerr}}_{k}|^{2}\approx|x^{\textit{in}}_{k}|^{2}-\gamma
L\operatorname{Im}\left(\sum_{\begin{subarray}{c}l-m+n=k\\\ n\neq k;l\neq
k\end{subarray}}x_{l}^{\textit{in}}\overline{x_{m}^{\textit{in}}}x_{n}^{\textit{in}}\overline{x_{k}^{\textit{in}}}\right)\
.$ (10)

The average power in our system is $P_{0}\approx 0.5$ W which leads to
invoking only a weak Kerr nonlinearity. A numerical simulation of the comb
evolution as a function of input power and propagation length is reported in
Fig. 1. The figure shows that the comb line intensities are not linearly
varying over the relevant power range, therefore showing that the first order
expansion in $\gamma$ used in Eq. 10 does not describe adequately the
experimental situation, although it does describe adequately the initial
propagation. Moreover, the neglect of the dispersion is not justified
experimentally. Nevertheless, Eq. 10, and its generalizations given in App. A
provide a rough idea of how the Kerr effect mixes information in the fiber.

#### Phase Modulation for Information Mixing.

To build a bridge between our current and previous work, we tested whether a
phase modulation applied to the output of a fiber can further increase the
performance of our ELM. As described in Refs. [17, 18], the phase modulation
is performed at comb-line spacing frequency $\Omega.$ It induces frequency-
domain interference between the lines and, thus, constitutes a possible
mechanism for effective information mixing. Contrary to Kerr-nonlinearity
mixing, the mixing process by phase modulation is a linear process.

We define $x^{\textit{PM}}_{k}$ as the comb-line amplitudes at the output of a
phase modulator [17]:

$\sum_{k}x^{\textit{PM}}_{k}e^{-ik(\omega_{0}+k\Omega)t}=\sum_{k}x^{\textit{Kerr}}_{k}e^{-i(\omega_{0}+k\Omega)t}e^{-im\cos{(\Omega
t)}}\ .$ (11)

where the time-dependent phase $m\cos{(\Omega t)}$ is the result of the phase
modulation, with $m$ being the coefficient defining the intensity of the phase
modulation. Using the Jacobi-Anger expansion $e^{-im\cos{(\Omega
t)}}=\sum_{l}i^{l}J_{l}(m)e^{-i(k+l)\Omega t},$ we can rewrite modulated comb
amplitudes as Bessel functions:

$x^{\textit{PM}}_{k}=\sum_{l}x^{\textit{Kerr}}_{l}i^{k-l}J_{k-l}(m)$ (12)

with $J_{l}$ being $l$-th order Bessel function of the first kind. In the
current setup, we use $m\approx 2$.

##  3 Methods

###  3.1 Experimental Setup

Figure 2: Experimental setup of a fiber-based Extreme Learning Machine. CW:
continuous-wave laser; PM: phase modulator; EDFA: erbium-doped fiber
amplifier; PSF: programmable spectral filter SMF: single-mode fiber; VOA:
variable optical attenuator; PD: photodiode. The source of the radio frequency
(RF) signal driving PM1 and (optionally) PM2 is not represented.

Fig. 2 shows the experimental setup. The input frequency comb is generated by
phase modulating the output of a continuous-wave laser (CW) with a lithium-
niobate phase modulator (PM1). The CW wavelength of
$\lambda=1554.6\;\mathrm{nm}$ corresponds to comb’s central wavelength. PM1 is
driven by a radio-frequency (RF) signal at frequency $\Omega\approx
17\;\mathrm{GHz}$ and power $P1\approx 30$ dBm. $\Omega$ determines the comb
line spacing, while $P1$ determines the number of comb lines. The current
setup generates approximately 25 comb lines spanning a spectrum of $\sim 3$
nm. Note that at the high RF power used to drive PM1, the transformations Eqs.
(11, 12) do not describe exactly the effect of PM1, probably because of
distortions in the sinusoidal RF signal (second-order corrections could be
applied to account for this effect, see [17] for details).

We encode the input data using a programmable spectral filter (PSF1, Coherent
II-VI Waveshaper). A PSF wavelength-multiplexes the comb lines, applies
variable attenuation and/or phase shifts to them, and recombines them to a
single output comb. Using PSF1, we encode the input data by applying either
attenuation or a phase delay to each comb line (Sec. 3.2).

A series of two erbium-doped-fiber amplifiers (EDFA, Keopsys models CEFA-C-PB-
LP and CEFA-C-BO-HP) increase the input comb power up to $500$ mW (starting
from a power of approximately $10$ mW). The characteristics of the erbium-
doped fibers are unknown, as the design of the employed devices is not
publicly available. Using openly available data sheets of other EDFA
manufactures (Corning ER) and assuming that the characteristics of the built-
in amplifying fibers would be similar, we estimate the following parameters
for amplifiers we use: the total amplifying fiber length is $20$ meters, the
dispersion parameter $@1550$ nm is $\approx 18$ ps/nm$\cdot$km corresponding
to $\beta_{2}\approx-23$ ps2/km, whereas the nonlinear coefficient $\gamma$ is
$\approx 6.3$ W-1 km${}^{-1}.$ Accordingly, the dispersion-induced phase shift
between two neighboring comb lines is
$\left(\beta_{2}\Omega^{2}L/2\right)_{\textit{EDFA}}\approx 0.3$ rad, whereas
the nonlinear phase shift is
$\left(\gamma\int_{0}^{L}P(z)dz\right)_{\textit{EDFA}}\approx 0.015$ rad
(assuming a constant amplification distributed along the propagation length).
Optionally, the comb can also be injected in $L=540$ m of a standard telecom
single-mode fiber (SMF, Throlabs SMF-28-J9). The relevant characteristics of
the fiber are: attenuation $\approx 0.18$ dB/km, dispersion parameter $@1550$
nm $\approx 18$ ps/nm$\cdot$km corresponding to $\beta_{2}\approx-23$ ps2/km,
nonlinear coefficient $\gamma\approx 1.2$ W-1 km-1. With these parameters, the
dispersion-induced phase shift between two neighboring comb lines is
$\left(\beta_{2}\Omega^{2}L/2\right)_{\textit{SMF}}\approx 9.0$ rad, while the
nonlinear phase shift is $\left(\gamma PL\right)_{\textit{SMF}}\approx 0.3$
rad (assuming $0.5$ W of power at the input of the fiber).

A variable optical attenuator (VOA) placed after the fiber reduces the optical
power down to a safe value compatible with the rest of the circuit.
Optionally, a second phase modulator (PM2) can be used to introduce an
additional transformation of the comb leading to further information mixing.
PM2 is driven by the same radio-frequency signal as PM1. However, its output
power is $P2\approx 20$ dBm due to a different electric amplification.

The output circuit is composed of a second programmable spectral filter (PSF2,
Coherent II-VI Waveshaper) and a photodiode (PD). PSF2 can be used either to
select individual comb lines, by applying band-pass filters, or to apply an
attenuation mask to the comb. In the first case, PD measures the intensity of
the selected comb line. In the second case, PD measures the sum of the comb
line intensities, each one weighted by the attenuation set on the mask. This
operation mode is employed to apply output weights optically and is described
in Sec. 3.2.

From the ANN perspective, PSF1 encodes the input layer of the ELM. EDFA, SMF
(optionally) and PM2 (optionally) provide a mix of input information, thus,
generating the hidden layer. PSF2 and PD constitute the output layer.

###  3.2 Working procedure

#### Input Layer.

In software implementations of ELMs, the dimension of the input layer
$N_{\textit{in}}$ depends on the task and is, in general, equal to the number
of features $N_{\textit{feat}}$ of the processed dataset, i.e.
$N_{\textit{in}}=N_{\textit{feat}}$. In our scheme, the number of input layer
neurons coincides with the number of available comb lines. This number is
governed by the parameters of the experimental setup. However, it can be
bigger than the number of data features, i.e.
$N_{\textit{in}}>N_{\textit{feat}}.$ We deploy a comb with
$N_{\textit{in}}\approx 25$ initial lines, while the number of features
depends on the benchmark dataset. Here, we use datasets with feature numbers
of 4 or 15, $N_{\textit{feat}}\in[4,\,15],$ to test our ELM.  It is therefore
necessary to establish a mapping between $N_{\textit{feat}}$ features to be
processed and the $N_{\textit{in}}$ available input neurons. To do so, we
stretch the input-feature vector $\mathbf{x}$ by an integer $d\geq 1$ that is
optimized for every task ([17]). The transformation reads as:

$\mathbf{x}=(x_{1},x_{2},\ldots,x_{N_{\textit{feat}}-1},x_{N_{\textit{feat}}})\rightarrow(\underbrace{x_{1},x_{1},\ldots,x_{1}}_{\text{$d$
times}},x_{2},\ldots,x_{{N_{\textit{feat}}}-1},\,\underbrace{x_{N_{\textit{feat}}},x_{N_{\textit{feat}}},\ldots,x_{N_{\textit{feat}}}}_{\text{$d$
times}}).$ (13)

The dimension of the stretched vector can still be smaller than the number of
initially available comb lines. The remaining, i.e. unmodulated,
$N_{\textit{in}}-N_{\textit{feat}}\cdot d$ lines are left unmodified and are
free to co-propagate through the ELM. This approach proves beneficial as
compared to full attenuation of the unused comb lines because it allows for
more optical power to propagate in the system, which, in turn, increases the
data-transforming impact of the Kerr nonlinearity. Moreover, the initially
unmodulated lines acquire information during the propagation, which
contributes to the increase of the dimensionality in the hidden layer. Thus, a
trade-off is required: high values of $d$ result in a lower number of
unmodulated lines which can contribute to the reduction of the total optical
power. On the other hand, low values of $d$ result in a lower number of
modulated lines, which can reduce the system’s response to the input. For each
dataset used here, we found an optimal value of $d$.

We tested the encoding of data features in comb line intensities and in
phases. When encoding in intensities, the features were converted to
attenuations in the range of [-30 dB, 0 dB]. In the case of phase encoding,
the features were converted to phase delays in the range of [0 rad, $2\pi$
rad]. Both attenuation and phase delays are set by PSF1.

#### Output Layer.

The hidden layer $\mathbf{h}$ is given by the intensities of the comb lines at
the output of VOA (or at the output of PM2 when it is used). The intensity
measurements are realized by using PSF2 as a tunable band-pass filter that
selects one line at a time to be subsequently measured by PD. Once
$\mathbf{h}$ is recorded, a ridge regression model is employed to retrieve the
optimal output weights $\mathbf{W}_{\textit{out}}$.

The system allows for two weighting methods that we call “digital weighing”
and “optical weighting”. In the digital weighting mode, the multiplication
$\mathbf{W}_{\textit{out}}\cdot\mathbf{h}$ is performed on a traditional
computer. This requires that the full vector of hidden neuron values
$\mathbf{h}$ is measured for each inference. In the optical weighting mode,
PSF2 applies attenuation to all comb lines simultaneously, which results in PD
measuring the sum of the weighted comb line intensities. To apply both
positive and negative weights, two consecutive intensity measurements are
performed: the first applying only the positive weights, and the second
applying only the (absolute values of the) negative weights. We define $y^{+}$
and $y^{-}$ as intensity measurements obtained with positive and negative
weights, respectively. The output signal $y$ is generated as their linear
combination:

$y=C^{+}\cdot y^{+}+C^{-}\cdot y^{-}+C^{0},$ (14)

where $(C^{+},\,C^{-},\,C^{0})$ is a set of coefficients optimized for each
task. The operation in Eq. (14) is currently executed on a digital computer,
but it can be implemented fully in the analog domain. This same scheme for
optical weighting has been presented in [17] and [20]. We also reported the
design for a fully integrated optical output layer based on a similar
principle in Refs. [28, 29].

###  3.3 Experimental Evaluation.

The evaluation of the system performance has been conducted as follows. First,
we tested our photonic ELM on three publicly available classification
benchmark tasks: the Heart disease[30], the Wine[31], and the Iris[32]
classification task. In this phase, we compared many possible operating
conditions of the system (i.e. encoding the input in comb-line intensities vs.
comb line phases or mixing the information through Kerr effect with and
without additional phase modulation).

After having identified the best-performing configuration using classification
datasets, we transfer the system parameters to test its performance on a
regression task. The selected benchmark predicts real-estate values, based on
the dataset first presented in Ref. [33]. For every tested configuration, we
evaluated the system performance for both, digital and optical weighting. For
both types of weighting, we tested the system with and without an optional
propagation in a 540-meter-long SMF.

We also compared the performance of the photonic ELM with a software ELM
counterpart with 25 hidden neurons and a quadratic or quartic nonlinearity.
The quadratic nonlinearity is inspired by our previous work [17], where the
only existing nonlinearity is due to the output-layer photodiode. The quartic
nonlinearity is inspired by our model for fiber-output comb-line intensities
(Eq. (10) and App. A). In the case of the Wine and Iris datasets, we compared
the performance of the current implementation with the results experimentally
obtained in the previously reported scheme that uses only phase modulation for
information mixing [17].

For every benchmark, 70% of each dataset were used for training the output
weights and 30% for testing the predictions. The classification performance
was evaluated in terms of accuracy, i.e. the percentage of correctly
classified examples over their total number. The performance on the regression
task was evaluated by means of the root-mean-square error (RMSE):

$\textit{RMSE}=\sqrt{\frac{1}{N}\sum_{1}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}},$
(15)

where $N$ is the number of regressions made, $y_{i}$ is the predicted value
for the $i$-th regression and $\hat{y}_{i}$ is the expected value for the
$i$-th regression. For each benchmark test, we repeated the training procedure
100 times selecting different partitions in train and test datasets. We
recorded both, the average and standard deviation of the accuracy values.

##  4 Results

In this section, we report the experimental results measured over multiple
configurations of the system, as described in Sec. 3.3. To distinguish among
the possible configurations, we employ the expressions “amplitude” and “phase”
to specify the encoding mechanism (whether data is encoded in comb line
amplitudes or phases, respectively). Also, we use the expressions “EDFA only”,
“EDFA+SMF” and “EDFA+SMF+PM” to specify the configuration of the setup (cf.
Fig. 2).

The first set of measurements is focussed on three classification benchmarks
and serves the purpose of identifying the best encoding method (amplitude vs.
phase) as well as the best experimental configuration (“EDFA+SMF” vs.
“EDFA+SMF+PM”). Encoding in amplitude in the “EDFA+SMF’ configuration clearly
emerged as the best performing option (we discuss this in the next
paragraphs). For ease of readability, the results of this first set of
measurements are reported in App. B, Tab. 1.

The second set of measurements was executed to test in depth the system in the
best configuration (encoding in amplitude, “EDFA+SMF”). Initially, we thought
that the EDFA did not play any role in the information mixing process, but
only helped boost the power so that the actual information mixing could take
place in the 540-long SMF. However, a careful analysis of our data showed that
significant information mixing took place in the EDFA itself. We, thus,
extended the study to include the “EDFA only’ configuration. This second
measurement set included both, the classification and the regression
benchmarks. The achieved results are reported in Fig. 3. In this figure, we
also compare the results with what was obtained from the previous experiment
(in which the information mixing was performed by a PM only [17]). We also
compare them with the results from an equivalent software-based ELMs (with 25
hidden neurons and a quadratic or quartic nonlinearity). The results are
summarized in App. B, Tabs. 2 and 3.

Figure 3: Experimental results obtained by employing digital (blue) and
optical (red) weighting, encoding the input in amplitudes, and using
"EDFA+SMF" and "EDFA only" configurations. They represent average scores
achieved by randomly selecting the train and test subsets 100 times. The error
bars represent the standard deviation of the accuracy distribution.
Experimental results are compared to fairly equivalent software ELMs (25
hidden nodes, quadratic or quartic nonlinearity) and to a previous experiment
based on information mixing via phase modulation only [17].

The recorded performance depends on both, the input-data encoding scheme
(amplitude vs. phase) and the mixing technique (i.e. the experimental
configuration). The encoding scheme and mixing technique can be considered as
hyperparameters of the network and, as usual in machine learning, the best-
performing hyperparameters depend on the task. Nevertheless, we can derive two
general observations from the results. First, the encoding in amplitude
generally performs better than encoding in phase. This result supports our
model that shows that the information mixing works via the comb-line
intensities rather than their phases and is driven by FWM (Eq. 10). Second, an
additional phase modulation by PM2 is never useful and in some cases even
detrimental. Phase modulation generates only linear mixing between the neurons
(Eq. (12)). Apparently, the nonlinear propagation in the fiber already
introduces sufficient data mixing, thus additional linear mixing provided by
the phase modulation does not improve the performance. Apart from that, the
electro-optical phase modulator has a polarization-dependent loss (nominally 5
dB) that might introduce fluctuations of the output power since thermal
effects can rotate the polarization state. These fluctuations, when happening
during the training phase, might result in a (further) reduction of the
performance.

In digital weighting mode, the photonic ELM achieved accuracies of $83.7\%\pm
4.4\%$, $99.9\%\pm 0.5\%$, and $100\%\pm 0\%$ on the Heart diseases, Wine, and
Iris classification benchmarks respectively. In optical weighting mode, the
photonic ELM achieved accuracies of $89.6\%\pm 8.5\%$, $100\%\pm 0\%$, and
$100\%\pm 0\%$ on the Heart diseases, Wine, and Iris classification benchmarks
respectively. Concerning the regression, the photonic ELM achieved an RMSE of
$8.52\pm 1.65$ and $8.18\pm 2.06$ in digital and optical weighting modes
respectively. These results refer to the best-performing configuration for
each task, on the test fraction of the dataset.

In the best-performing configuration, optical weighing always delivers results
comparable to digital weighting. In some cases, optical weighting average
accuracies are slightly higher than digital-weighting ones. This is probably
due to some statistical fluctuations. However, we note that, in optical-
weighting mode, the positive-weighted signals and the negative-weighted
signals are combined through optimized coefficients (Eq. (14)). In principle,
this combination should be linear. In practice, though, an extra nonlinearity
may be introduced by the saturation of the photodiodes, which might explain
the increased performance.

When compared to our previous experiment (Fig. 3 and App. B, Tab. 3), the
approach based on the amplitude-encoded data mixed by the Kerr nonlinearity
delivers better performance in both (digital and optical) weighting modes.

Our experimental scheme also slightly outperformed the software-based ELM
counterparts that were designed to replicate the experimental information
mixing (with quadratic or quartic nonlinearity). This indicates that the
experimental ELM might benefit from nonlinearities that are stronger and/or
richer than the quadratic and quartic ones. This could be due to higher-order
contributions of the Kerr effect that are not captured by the quartic
approximation, as well as nonlinearities in other experimental devices (e.g.,
the saturation of the photodiodes).

Further, we observed that propagation only in the erbium-doped fiber is enough
to perform information processing, although the addition of the 540-meter-long
SMF systematically results in a slightly increased performance. To better
clarify the effect of the propagation, we numerically integrate Eq. (7) to
evaluate both the performance and accumulated nonlinear phase as a function of
propagation length. The results are reported in Fig. 4. They show that the
propagation in the EDFA ($L<20$ m) is sufficient for the ELM performance to
approach the best value. We reduced the plots to a length of $40$ meters
because this is the most interesting region, no performance increase was
detected for longer propagation distances.

Figure 4: Results of the numerical simulation based on Eq. (7) for two
benchmark tasks: Iris classification (a) and Real Estate Value regression (b).
The initial $20$ meters of propagation take place in an erbium-doped
(amplifying) fiber, while the last $20$ meters happen in a standard telecom
fiber (Sec. 3.1). The cumulated nonlinear phase (left axis) is given by
$\Phi(L)=\int_{0}^{L}\gamma(z)P(z)dz$ where the $z$ dependence of $\gamma$
takes into account that the nonlinear coefficient is different in the EDFA and
in the standard fiber, while the $z$ dependence of $P(z)$ takes into account
amplification and loss. The discontinuity of $\Phi(L)$ at $L=20$ m is because
the erbium-doped fiber is amplifying, thus, $P(z)$ exponentially increases
with length, while the standard telecom fiber is not amplifying, thus, $P(z)$
slightly decreases due to attenuation. The results of this figure were
obtained using the 4th-order Runge-Kutta in the Interaction Picture Method
(RK4IP) including a realistic description of the erbium-doped fiber [34].

##  5 Conclusion

We presented a photonic implementation of a frequency-multiplexed photonic
Extreme Learning Machine that encodes input information in the lines of a
frequency comb and exploits the Kerr nonlinearity in optical fibers to process
it. In our system, the information processing occurs via Kerr-driven four-wave
mixing, rather than via self- or cross-phase modulation or dispersion effects,
which is shown by a perturbative model we derived and presented here. We
experimentally validated the Extreme Learning Machine using four different
benchmark tasks including classification and regression problems. Optimizing
the system performance, we compared multiple possible configurations, varying
the input encoding method (amplitude vs. phase), way of information mixing
(Kerr nonlinearity with and without a subsequent phase modulation), and the
kind of output weights (digital vs. optical). The best results were achieved
by encoding the neurons in amplitudes of the frequency-comb lines and using
only the Kerr for information mixing. Digital and optical weighting provided
comparable results. Encoding the neurons in the phases of the comb lines
resulted in a worsening of the EML performance.

Surprisingly, our work shows that a weak Kerr nonlinearity can yield good
results on classification and regression tasks. In fact the propagation in the
$20$-meter-long erbium-doped fiber included in the amplification stage is
already sufficient to process information. The amount of nonlinearity can be
quantified by the product $\phi=\gamma LP$, which in our work is either
$\phi\approx 0.0015$ rad or $\phi\approx 0.3$ rad depending on the
configuration, corresponding to a very weak nonlinearity. All our results,
theoretical, experimental, numerical-simulations, coincide well with each
other showing that a weak Kerr nonlinearity is a highly efficient mechanism
for information mixing in our frequency-multiplexed fiber-based Extreme
Learning Machine.

Let us now briefly compare our Extreme Learning Machine with approaches
presented in Refs. [21] and [23]. In these references, the authors use broad
continuous spectra from femtosecond pulsed lasers as sources to modulate the
input data onto them by deploying programmable spectral filters. Then, they
propagate modulated spectra in pieces of highly nonlinear fibers to allow for
data transformation and collect fiber outputs for further evaluation on a PC.
Both papers point out the importance of the nonlinear dynamics described by
the Nonlinear Schrödinger equation. Intuitively, they expect that higher
values of optical Kerr nonlinearity would lead to a better performance. This
is why their setups operate in highly nonlinear regimes (for instance in Ref.
[21] the nonlinear phase is $\phi\approx 165$ rad). This intuition is
justified as theoretical studies from Ref. [8] suggest. Contrary to it, we
show that already a weak Kerr nonlinearity can be used for effective data
processing.

As the information-carrying spectra from Refs. [21] and [23] are continuos,
they are subject to several effects such as, for instance, dispersion, self-
phase modulation, intrapulse stimulated Raman scattering, and four-wave
mixing. All these effects occur simultaneously and impact the dynamics in a
complex, delicate, and interconnected manner. In our frequency-multiplexed
Extreme Learning Machine operating with a discrete spectrum under weak Kerr
nonlinearity, the only really decisive effect is four-wave mixing, while
dispersion plays a minor role. In general, the dynamics of systems that
operate in highly nonlinear regime are sensitive to the input noise. This
dependence can severely affect the data-processing capability of a
neuromorphic scheme. From this perspective, weak-nonlinearity schemes like
ours can offer the possibility of more robust neuromorphic-computing
implementations. One way or another, both approaches (with strong or weak Kerr
nonlinearity) have advantages and disadvantages and should be pursued in
future. For both approaches, it is important to understand the evolution of
the data-encoded optical signal and how various effects might enhance or
impact the overall data-processing capability of the neuromorphic scheme.

In future, more complex task with more input features should be addressed. For
this, we would need a broader frequency comb to encode the data. A broader
comb can be generated using the beating between two continuous-wave lasers,
rather than deploying a phase modulator at the initial stage (cf. [35]). To
increase the data throughput, we will replace the quite slow programmable
spectral filters we currently deploy by photonic-integration designs we
present in Refs. [28], [29]. In fact, our whole Extreme Learning Machine can
be transferred to a chip. Semiconductor-based solutions needed for on-chip
integration are easily available (e.g. semiconductor optical amplifiers incl.
erbium-doped waveguides [36] as well as frequency-comb sources [37], [38]).
Photonic integration of our Extreme Learning Machine is attractive as it will
allow for lowering the power consumption while contributing to scalability,
stability, and practicability.

## Acknowledgement

We acknowledge interesting and inspiring conversations with Sonia Boscolo
(Aston University Birmingham, UK), Christophe Finot (ICB - Université de
Bourgogne, France), and Ingo Fischer (IFISC - Universitat de les Illes
Balears, Spain).

## Funding

We acknowledge the funding from H2020 Marie Skłodowska-Curie Actions (Project
POST-DIGITAL, grant number 860830), FWO and F.R.S.-FNRS Excellence of Science
(EOS) program grant 40007536, and Win4Space project of Walloon’s Region
excellence programme.

## Disclosures

The authors declare no conflicts of interest.

##  Appendix A Perturbative approximation of the Kerr effect on comb line
intensities

###  A.1 Nonlinear Schrödinger Equation

The propagation along the fiber is described by the Nonlinear Schrödinger
Equation (NLS) [27]:

$\partial_{z}x(z,t)=i\beta(i\partial_{t})x(z,t)-\frac{\alpha}{2}x(z,t)-i\gamma|x(z,t)|^{2}x(z,t),$
(16)

where $x(z,t)$ is the electric field, $\beta(\omega)$ the propagation
constant, $\alpha$ the attenuation (negative $\alpha$ corresponds to the
gain), and $\gamma$ the nonlinear Kerr coefficient of the fiber.

The total power at the propagation point $z$ is given by $P(z)=\int
dt|x(z,t)|^{2}.$ If $\beta(i\partial_{t})$ is a Hermitian operator, then the
total power satisfies the equation$\frac{dP(z)}{dz}=-\alpha P(z)$. Therefore,
we have:

$P(z)=P_{0}e^{-\alpha z}$ (17)

with $P_{0}$ being the initial total power.

According to Sec. 2, the input field, thus, the initial condition, is a
frequency comb:

$x(z=0,t)=\sum_{k}x_{k}e^{-i\Omega kt},$ (18)

where $x_{k}$ are the (complex) amplitudes of the input comb lines and
$\Omega$ is the spacing between the comb lines.

If $\gamma=0$, then the solution of the NLS reads as:

$x(z,t)=\sqrt{P_{0}}\sum_{k}b_{k}^{0}e^{-i\Omega kt}e^{i\beta(\Omega
k)z}e^{-\frac{\alpha}{2}z},$ (19)

which shows that each comb line acquires a different phase due to the
dispersion ($\beta(\Omega k)z$) and all comb lines are attenuated identically
($-\frac{\alpha}{2}z$).

If $\gamma\neq 0$, we write the solution as:

$\displaystyle x(z,t)$ $\displaystyle=$
$\displaystyle\sum_{k}x_{k}(z)e^{-i\Omega kt}$ (20) $\displaystyle=$
$\displaystyle\sqrt{P_{0}}e^{-\frac{\alpha}{2}z}\sum_{k}b_{k}(z)e^{-i\Omega
kt}e^{i\beta(\Omega k)z},$

which shows that the temporal comb shape is preserved. However, the amplitudes
and phases of the comb lines acquire a $z-$dependence due to the Kerr
nonlinearity. Here, $b_{k}(z)$ are normalized comb-line amplitudes, i.e.

$\sum_{k}|b_{k}(z)|^{2}=1,\ $ (21)

to satisfy Eq. (17).

Inserting Eq. (20) in Eq. (16), we obtain:

$\sum_{k}\partial_{z}b_{k}e^{-i\Omega kt}e^{i\beta(\Omega
k)z}e^{-\frac{\alpha}{2}z}=-i\gamma
P_{0}\sum_{l,m,n}b_{l}\overline{b}_{m}b_{n}e^{-i\Omega(l-m+n)t}e^{i\left[\beta(\Omega
l)-\beta(\Omega m)+\beta(\Omega n)\right]z}e^{-\frac{3\alpha}{2}z}.$ (22)

We can rewrite Eq. 22 as equations for each comb line amplitude:

$\partial_{z}b_{k}=-i\gamma P_{0}e^{-\alpha
z}\sum_{l-m+n=k}b_{l}\overline{b}_{m}b_{n}e^{i\left[\beta(\Omega
l)-\beta(\Omega m)+\beta(\Omega n)-\beta(\Omega k)\right]z}\ .$ (23)

Fig. 1(a) shows the comparison between the numerical integration of Eq. (23)
and the NLS integration via the split-step Fourier method (SSFM). The results
deviate only slightly, which is probably due to a limited numerical precision.
Thus, we conclude that Eq. 23 represents a good description of the comb-line
evolution in optical fibers.

###  A.2 Self- and Cross-Phase Modulation

The $z-$dependent phase on the right-hand side of Eq. (23) vanishes if the sum
indices $l$ or $n$ equal $k$. These terms correspond to self- and cross-phase
modulation. If we restrict the right-hand side of Eq. (23) to the terms such
that $l$ or $n$ equal $k,$ we have the equations:

$\displaystyle\partial_{z}b_{k}$ $\displaystyle=$ $\displaystyle-i\gamma
P_{0}e^{-\alpha z}\left(2\sum_{l\neq
k}|b_{l}|^{2}+|b_{k}|^{2}\right)b_{k}e^{-\alpha z}$ (24) $\displaystyle=$
$\displaystyle-i\gamma P_{0}e^{-\alpha z}(2-|b_{k}|^{2})b_{k},$

where we used Eq. (21). Eq. (24) is solved by:

$b_{k}(z)=e^{-i\phi_{k}(z)}b_{k}^{0},$ (25)

with

$\phi_{k}(z)=\gamma P_{0}(2-|b_{k}^{0}|^{2})\frac{1-e^{-\alpha z}}{\alpha}\ .$
(26)

Here, $b_{k}^{0}$ is the initial normalized comb-line amplitude at $z=0$. Eq.
(25) shows that only the phase of $b_{k}$ is affected by self- and cross-phase
modulation, not its amplitude. In our implementation of a photonic ELM, self-
and cross-phase modulation have no measurable effect as we read out only the
comb intensities.

###  A.3 Four-Wave Mixing

Eq. (24) neglects FWM terms in Eq. (23). To take them into account, we carry
out a perturbative expansion in $\gamma$. For this, we perturb the comb-line
amplitudes by a small function $\epsilon_{k}(z):$

$x(z,t)=\sqrt{P_{0}}\sum_{k}e^{-i\Omega kt}e^{i\beta(\Omega
k)z}e^{-i\phi_{k}(z)}e^{-\frac{\alpha}{2}z}\left(b_{k}^{0}+\gamma\epsilon_{k}(z)\right),$
(27)

where $\phi_{k}(z)$ given by Eq. (26) and takes into account the self- and
cross-phase modulation.

Upon inserting Eq. (27) into the NLS (Eq. (16)), the terms of order
$\gamma^{0}$ cancel, while the terms of order $\gamma^{1}$ give an equation
for $\epsilon_{k}$(z):

$\partial_{z}\epsilon_{k}=-iP_{0}e^{-\alpha z}\sum_{l-m+n=k;n\neq k;l\neq
k}b_{l}^{0}\overline{b_{m}^{0}}b_{n}^{0}e^{i\Delta\beta_{k,l,m,n}z}+O(\gamma)$
(28)

where the condition on the summation indices $l$, $m$, and $n$ expresses that
self- and cross-phase terms are excluded, and we denote

$\Delta\beta_{k,l,m,n}=\beta(\Omega l)-\beta(\Omega m)+\beta(\Omega
n)-\beta(\Omega k)\ .$ (29)

Note that we have omitted the contribution of the phases $\phi_{k}(z)$ in Eq.
(28) as these are of order $\gamma$.

The solution of Eq. (28) gives us the perturbation function $\epsilon_{k}(z):$

$\epsilon_{k}(z)=P_{0}\sum_{l-m+n=k;n\neq k;l\neq
k}\frac{b_{l}^{0}\overline{b_{m}^{0}}b_{n}^{0}}{\Delta\beta_{k,l,m,n}+i\alpha}\left(1-e^{i\Delta\beta_{k,l,m,n}z}e^{-\alpha
z}\right).$ (30)

It will allow us to determine the intensities of the comb lines in the next
step.

###  A.4 Comb-line Intensities

The intensity of a comb-line $k$ is obtained by taking the norm square of Eq.
(27) to yield:

$I_{k}=|x_{k}(z)|^{2}=P_{0}e^{-\alpha z}\left(|b_{k}^{0}|^{2}\\\
+\gamma(\overline{b_{k}^{0}}\epsilon_{k}(z)+c.c.)+O(\gamma^{2})\right)\ .$

The change in the intensity of this line to order $\gamma^{1}$ (noted $\Delta
I_{k}^{1}$), normalized by the total power $P(z)=P_{0}e^{-\alpha z}$ is

$\displaystyle\frac{\Delta I_{k}^{1}}{P(z)}$ $\displaystyle=$
$\displaystyle\gamma(\overline{b_{k}^{0}}\epsilon_{k}+c.c.)$ (31)
$\displaystyle=$ $\displaystyle\gamma P_{0}\ Real\left[\sum_{l-m+n=k;n\neq
k;l\neq
k}\frac{b_{l}^{0}\overline{b_{m}^{0}}b_{n}^{0}\overline{b_{k}^{0}}}{\Delta\beta_{k,l,m,n}+i\alpha}\left(1-e^{i\Delta\beta_{k,l,m,n}z}e^{-\alpha
z}\right)\right].$

A comparison between the linearized change in intensity (Eq. (31)) and the
result of the numerical integration of NLS is given in Fig. 1(b). With
increasing propagation length, we see a deviation between the solution
provided by NLS and Eq. (31).

Eq. (31) simplifies in the limit of small propagation lengths when
$e^{i\Delta\beta_{k,l,m,n}z}e^{-\alpha z}\approx
1+(i\Delta\beta_{k,l,m,n}-\alpha)z,$ whereupon one finds:

$\frac{\Delta I_{k}^{1}}{P(z)}\approx-\gamma P_{0}z\
Im\left[\sum_{l-m+n=k;n\neq k;l\neq
k}b_{l}^{0}\overline{b_{m}^{0}}b_{n}^{0}\overline{b_{k}^{0}}\right].$ (32)

Here, we see that the nonlinearity that is introduced to the ELM network by
the Kerr effect is of the fourth order, i.e. quartic.

Eq. (31) also simplifies in the case of strong amplification when $e^{-\alpha
z}\gg 1$ (corresponding to propagation in the EDFA). Assuming also that
$\alpha\gg\Delta\beta$, one finds

$\frac{\Delta I_{k}^{1}}{P(z)}\approx\gamma P_{0}e^{-\alpha z}\
Im\left[\sum_{l-m+n=k;n\neq k;l\neq
k}b_{l}^{0}\overline{b_{m}^{0}}b_{n}^{0}\overline{b_{k}^{0}}e^{i\Delta\beta_{k,l,m,n}z}\right].$
(33)

##  Appendix B Summary of the Results

Here we present all our results in three tables.

|  | Amplitude Encoding | Phase Encoding
---|---|---|---
Task |  | Weighting
---
EDFA+SMF |  | EDFA+SMF+PM
---
EDFA+SMF |  | EDFA+SMF+PM
---
Heart disease | Digital | $83.7\%$ | $80.5\%$ | $71.5\%$ | $60.6\%$
Optical | $89.6\%$ | $83.3\%$ | $80.4\%$ | $56.1\%$
Wine | Digital | $99.9\%$ | $99.3\%$ | $92.8\%$ | $96.7\%$
Optical | $100\%$ | $91.5\%$ | $95.4\%$ | $78.0\%$
Iris | Digital | $100\%$ | $100\%$ | $100\%$ | $99.8\%$
Optical | $100\%$ | $99.8\%$ | $78.2\%$ | $87.7\%$
Table 1: Results for the first measurement set. Three classification tasks were used in four experimental schemes. From these results, it appears that the amplitude encoding of input information is better than the phase encoding. Also, we see that using a phase modulator for information mixing does not improve the performance of the ELM. Task | Weighting | EDFA+SMF | EDFA only
---|---|---|---
(Amplitude encoding) | (Amplitude encoding)
Heart disease | Digital | 83.7% | 82.6%
(Accuracy) | Optical | 89.6% | 83.2%
Wine | Digital | 99.9% | 100%
(Accuracy) | Optical | 100% | 99.6%
Iris | Digital | 100% | 98.6%
(Accuracy) | Optical | 100% | 88.6%
Real estate | Digital | 8.52 | 8.91
(RMSE) | Optical | 8.96 | 8.18
Table 2: Results for the second measurement set. All benchmarks were executed for the best performing configuration selected after the first measurement set (amplitude encoding, without phase modulation). Each task was evaluated with and without the propagation in a 540-meter-long single-mode fiber. Adding an SMF after the EDFA generally results in a small performance improvement. Task | Software ELM | Previous experiment
---|---|---
Quadratic NL | Quartic NL | (Amplitude, PM only)
| Heart disease
---
(Accuracy)
81.9% | 79.6% | Not measured
| Wine
---
(Accuracy)
97.5% | 96.2% | 97.5%
| Iris
---
(Accuracy)
97.3% | 93.5% | 93.9%
| Real estate
---
(RMSE)
8.79 | 8.73 | Not measured
Table 3: Results obtained with a software ELM as well as experimental results
from the previous work where information mixing was facilitated by a phase
modulator [17].

## References

  * [1] “Beyond von Neumann,” Nature Nanotechnology 15, 507–507 (2020).
  * [2] A. Reuther, P. Michaleas, M. Jones, V. Gadepally, S. Samsi, and J. Kepner, “Survey and benchmarking of machine learning accelerators,” in _2019 IEEE High Performance Extreme Computing Conference (HPEC),_ (IEEE, 2019).
  * [3] D. Marković, A. Mizrahi, D. Querlioz, and J. Grollier, “Physics for neuromorphic computing,” Nature Reviews Physics 2, 499–510 (2020).
  * [4] K. Nakajima, “Physical reservoir computing—an introductory perspective,” Japanese Journal of Applied Physics 59, 060501 (2020).
  * [5] G.-B. Huang, Q.-Y. Zhu, and C.-K. Siew, “Extreme learning machine: theory and applications,” Neurocomputing 70, 489–501 (2006).
  * [6] G.-B. Huang, “What are extreme learning machines? filling the gap between frank rosenblatt’s dream and john von neumann’s puzzle,” Cognit. Comput. 7, 263–278 (2015).
  * [7] G.-B. Huang, L. Chen, and C.-K. Siew, “Universal approximation using incremental constructive feedforward networks with random hidden nodes,” IEEE Transactions on Neural Networks 17, 879–892 (2006).
  * [8] G. Marcucci, D. Pierangeli, and C. Conti, “Theory of neuromorphic computing by waves: Machine learning by rogue waves, dispersive shocks, and solitons,” Physical Review Letters 125 (2020).
  * [9] B. J. Shastri, A. N. Tait, T. F. de Lima, W. H. P. Pernice, H. Bhaskaran, C. D. Wright, and P. R. Prucnal, “Photonics for artificial intelligence and neuromorphic computing,” Nature Photonics 15, 102–114 (2021).
  * [10] A. Saade, F. Caltagirone, I. Carron, L. Daudet, A. Dremeau, S. Gigan, and F. Krzakala, “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),_ (IEEE, 2016).
  * [11] D. Pierangeli, G. Marcucci, and C. Conti, “Photonic extreme learning machine by free-space optical propagation,” Photonics Res. 9, 1446–1454 (2021).
  * [12] U. Teğin, M. Yıldırım, İ. Oğuz, C. Moser, and D. Psaltis, “Scalable optical learning operator,” Nature Computational Science 1, 542–549 (2021).
  * [13] S. Ortín, M. C. Soriano, L. Pesquera, D. Brunner, D. San-Martín, I. Fischer, C. Mirasso, and J. Gutiérrez, “A unified framework for reservoir computing and extreme learning machines based on a single time-delayed neuron,” Sci. Rep. 5, 14945 (2015).
  * [14] X. Xu, M. Tan, B. Corcoran, J. Wu, T. G. Nguyen, A. Boes, S. T. Chu, B. E. Little, R. Morandotti, A. Mitchell, D. G. Hicks, and D. J. Moss, “Photonic perceptron based on a kerr microcomb for high-speed, scalable, optical neural networks,” Laser Photonics Rev. 14, 2000070 (2020).
  * [15] X. Xu, M. Tan, B. Corcoran, J. Wu, A. Boes, T. G. Nguyen, S. T. Chu, B. E. Little, D. G. Hicks, R. Morandotti, A. Mitchell, and D. J. Moss, “11 tops photonic convolutional accelerator for optical neural networks,” Nature 589, 44–51 (2021).
  * [16] J. Feldmann, N. Youngblood, M. Karpov, H. Gehring, X. Li, M. Stappers, M. L. Gallo, X. Fu, A. Lukashchuk, A. S. Raja, J. Liu, C. D. Wright, A. Sebastian, T. J. Kippenberg, W. H. P. Pernice, and H. Bhaskaran, “Parallel convolutional processing using an integrated photonic tensor core,” Nature 589, 52–58 (2021).
  * [17] A. Lupo, L. Butschek, and S. Massar, “Photonic extreme learning machine based on frequency multiplexing,” Optics Express 29, 28257 (2021).
  * [18] A. Lupo and S. Massar, “Parallel extreme learning machines based on frequency multiplexing,” Applied Sciences 12 (2022).
  * [19] M. Yildirim, I. Oguz, F. Kaufmann, M. R. Escalé, R. Grange, D. Psaltis, and C. Moser, “Nonlinear optical feature generator for machine learning,” APL Photonics 8 (2023).
  * [20] L. Butschek, A. Akrout, E. Dimitriadou, A. Lupo, M. Haelterman, and S. Massar, “Photonic reservoir computer based on frequency multiplexing,” Optics Letters 47, 782 (2022).
  * [21] T. Zhou, F. Scalzo, and B. Jalali, “Nonlinear schrödinger kernel for hardware acceleration of machine learning,” Journal of Lightwave Technology 40, 1308–1319 (2022).
  * [22] M. Sorokina, “Multi-channel optical neuromorphic processor for frequency-multiplexed signals,” Journal of Physics: Photonics 3, 014002 (2020).
  * [23] B. Fischer, M. Chemnitz, Y. Zhu, N. Perron, P. Roztocki, B. MacLellan, L. D. Lauro, A. Aadhi, C. Rimoldi, T. H. Falk, and R. Morandotti, “Neuromorphic computing via fission-based broadband frequency generation,” Advanced Science (2023).
  * [24] J. Pauwels, G. Verschaffelt, S. Massar, and G. Van der Sande, “Distributed kerr non-linearity in a coherent all-optical fiber-ring reservoir computer,” Frontiers in Physics 7 (2019).
  * [25] B. Kryzhanovsky, L. Litinskii, and A. Fonarev, “Parametrical neural network based on the four-wave mixing process,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 502, 517–519 (2003).
  * [26] B. Kryzhanovsky, V. Kryzhanovsky, A. Mikaelian, and A. Fonarev, “Parametric dynamic neural network recognition power,” Optical Memory and Neural Networks 10, 211–218 (2001).
  * [27] G. Agrawal, _Nonlinear Fiber Optics, Fifth ed._ (Elsevier, 2013).
  * [28] T. Jonuzi, A. Lupo, M. C. Soriano, S. Massar, and J. D. Domenéch, “Integrated programmable spectral filter for frequency-multiplexed neuromorphic computers,” Optics Express 31, 19255 (2023).
  * [29] T. Jonuzi, A. Lupo, M. C. Soriano, J. D. D. Gomez, and S. Massar, “Integrated optical output layer for a reservoir computer based on frequency multiplexing,” in _AI and Optical Data Sciences IV,_ K. ichi Kitayama and B. Jalali, eds. (SPIE, 2023).
  * [30] A. Janosi, W. Steinbrunn, M. Pfisterer, and R. Detrano, “Heart diseases data set,” https://archive.ics.uci.edu/ml/datasets/heart+disease.
  * [31] M. Forina, “Wine data set,” https://archive.ics.uci.edu/ml/datasets/wine.
  * [32] R. Fisher, “Iris data set,” http://archive.ics.uci.edu/ml/datasets/iris.
  * [33] I.-C. Yeh and T.-K. Hsu, “Building real estate valuation models with comparative approach through case-based reasoning,” Applied Soft Computing 65, 260–271 (2018).
  * [34] S. Balac, A. Fernandez, F. Mahé, F. Méhats, and R. Texier-Picard, “The interaction picture method for solving the generalized nonlinear schrödinger equation in optics,” ESAIM: Mathematical Modelling and Numerical Analysis 50, 945–964 (2016).
  * [35] M. Zajnulina, J. M. Chavez Boggio, M. Böhm, A. A. Rieznik, T. Fremberg, R. Haynes, and M. M. Roth, “Generation of optical frequency combs via four-wave mixing processes for low- and medium-resolution astronomy,” Applied Physics B 120, 171–184 (2015).
  * [36] Y. Liu, Z. Qiu, X. Ji, A. Lukashchuk, J. He, J. Riemensberger, M. Hafermann, R. N. Wang, J. Liu, C. Ronning, and T. J. Kippenberg, “A photonic integrated circuit–based erbium-doped amplifier,” Science 376, 1309–1313 (2022).
  * [37] M. Zajnulina, B. Lingnau, and K. Lüdge, “Four-wave mixing in quantum-dot semiconductor optical amplifiers: A detailed analysis of the nonlinear effects,” IEEE Journal of Selected Topics in Quantum Electronics 23, 1–12 (2017).
  * [38] A. L. Gaeta, M. Lipson, and T. J. Kippenberg, “Photonic-chip-based frequency combs,” Nature Photonics 13, 158–169 (2019).