# Machine learning for online sea ice bias correction within global ice-ocean
simulations

# Supporting Information: Machine learning for online sea ice bias correction
within global ice-ocean simulations

###### Abstract

In this study we perform online sea ice bias correction within a GFDL global
ice-ocean model. For this, we use a convolutional neural network (CNN) which
was developed in a previous study [Gregory . (2023)] for the purpose of
predicting sea ice concentration (SIC) data assimilation (DA) increments. An
initial implementation of the CNN shows systematic improvements in SIC biases
relative to the free-running model, however large summertime errors remain. We
show that these residual errors can be significantly improved with a data
augmentation approach, in which sequential CNN and DA corrections are applied
to a new simulation over the training period. This then provides a new
training data set with which to refine the weights of the initial network. We
propose that this machine-learned correction scheme could be utilized for
generating improved initial conditions, and also for real-time sea ice bias
correction within seasonal-to-subseasonal sea ice forecasts.

Geophysical Research Letters

Atmospheric and Oceanic Sciences Program, Princeton University, NJ, USA
Geophysical Fluid Dynamics Laboratory, NOAA, Princeton, NJ, USA Courant
Institute of Mathematical Sciences, New York University, New York, NY, USA

Will<EMAIL_ADDRESS>
We use a convolutional neural network (CNN) to perform online sea ice bias
correction within global ice-ocean simulations.

The CNN systematically reduces the free-running model bias in both the Arctic
and Antarctic.

The online performance can be improved by combining CNN and data assimilation
corrections in order to iteratively augment the training data.

## Plain Language Summary

Climate models contain errors which often lead to predictions which are
consistently out of agreement with what we observe in reality. In some cases
we know the origin of these errors, for example predicting too much sea ice as
a result of consistently cool ocean temperatures. In reality however, there
are typically numerous model errors interacting across the atmosphere, ocean
and sea ice, and to manually parse through large volumes of climate model data
in an attempt to isolate these errors in time and space is highly impractical.
Machine learning on the other hand is a framework which is well-suited to this
task. In this work we take a machine learning model which, at any given
moment, ingests information about a climate model’s atmosphere, ocean and sea
ice conditions, and predicts how much error there is in the climate model’s
representation of sea ice, without seeing any actual sea ice observations. We
use this to adjust the sea ice conditions in one particular climate model as
it is running forward in time making predictions, and we find that this
significantly reduces the model’s sea ice errors globally.

##  1 Introduction

Machine learning (ML) algorithms are beginning to cement their position as
viable subgrid-scale climate model parameterizations, through their ability to
isolate complex non-linear relationships within large volumes of high
dimensional data [Brenowitz Bretherton (2018), Gentine . (2018), O’Gorman
Dwyer (2018), Yuval O’Gorman (2020), Finn . (2023), Sane . (2023)]. Typically
this is achieved by training an ML model to learn a functional mapping which
characterizes the impact of subgrid processes on resolved scales, by training
on high resolution simulations or observational data. Significant effort is
currently being afforded to the development of these ML parameterizations in
the context of e.g., ocean turbulence, with early results [Zanna Bolton
(2020), Frezat . (2022), Ross . (2023), Kurz . (2023), C. Zhang . (2023)]
highlighting their potential to improve important climate statistics, such as
eddy kinetic energy at large scales, over their traditional physics-based
counterparts.

Alternatively, combining data assimilation (DA) and ML has shown to be a
promising framework for learning either subgrid parameterizations or
systematic model errors across various domains [Bonavita Laloyaux (2020),
Brajard . (2021), Farchi . (2021), Mojgani . (2022), Chen . (2022), Laloyaux .
(2022), He . (2023)]. In a recent study by Gregory2023, hereafter G23, the
authors presented a DA-based ML framework in which convolutional neural
networks (CNNs) were used to predict state-dependent sea ice errors within an
ice-ocean configuration of the Geophysical Fluid Dynamics Laboratory (GFDL)
Seamless system for Prediction and EArth System Research (SPEAR) model, as a
way to highlight the feasibility of a data-driven sea ice model
parameterization within SPEAR. They approached this by first showing that the
climatological sea ice concentration analysis increments from an ice-ocean DA
experiment map closely onto the systematic bias patterns of the equivalent
free-running model. This suggested that an ML model which is able to predict
the analysis increments could, in principle, reduce sea ice biases as an
online model parameterization or bias correction tool. Their subsequent CNN
architecture then used information from local model state variables and their
tendencies, to make predictions of the corresponding sea ice concentration
analysis increment at any grid cell location. These offline predictions were
shown to generalize well to both the Arctic and Antarctic domains, and across
all seasons. However, offline performance does not always directly translate
to online simulations, which can sometimes exhibit instabilities as well as
climate drift after implementation [Rasp . (2018), Ott . (2020), Brenowitz .
(2020)]. In such cases, the ML model may require an additional online training
step in order to sample a larger model state space to which it was initially
trained [Rasp (2020)].

In this present work, we advance the field of ML-based parameterizations by
investigating the online performance of the G23 DA-based ML model when used as
a tool to correct short-term sea ice error growth. We implement the correction
scheme here within a coupled ice-ocean configuration of SPEAR, as the G23 CNN
was originally trained on data from an ice-ocean DA system, which therefore
allows us to make direct comparisons of model biases and increments produced
from both the CNN and DA simulations. If the CNN is able to reduce sea ice
biases relative to the free-running model, then this will provide a solid
foundation for future work into assessing the generalization to fully coupled
systems, and ultimately a physics-based sea ice model parameterization.

##  2 Data and methods

###  2.1 SPEAR ice-ocean model

SPEAR is a fully coupled ice-ocean-atmosphere-land model [Delworth . (2020)],
which shares the same components as the GFDL CM4 model [Held . (2019)],
however with parameterizations and resolutions geared toward seasonal-to-
decadal prediction. The ocean and sea ice components are configured at a 1∘
horizontal resolution and correspond to the Modular Ocean Model v6 (MOM6) and
the Sea Ice Simulator v2 (SIS2), respectively [Adcroft . (2019)]. In this
work, we consider an ice-ocean configuration of SPEAR, in which MOM6 and SIS2
are forced by atmospheric conditions from the Japanese 55-year Reanalysis for
driving ocean-sea-ice models (JRA55-do; Tsujino2018). Details of the ice-ocean
experiments are provided in section 2.3.

Figure 1: Schematic of bias correction schemes, shown in each panel for one
5-day assimilation/correction cycle. The dots represent the daily model state
integrating forward in time. The vertical lines are then the corrections from
either the CNN or DA. (a) Out-of-the-box G23 CNN training and implementation
showing the Z21 DA simulation in red and the CNN simulation in light blue. (b)
Optimized G23 CNN training and implementation showing simulations with
combined CNN and DA corrections in violet, and the optimized CNN simulation in
dark blue.

###  2.2 Machine learning model

The CNN model from G23 was trained to predict sea ice concentration (SIC)
increments from a SPEAR ice-ocean DA experiment (Zhang2021; hereafter Z21).
The Z21 DA experiment spanned January 1st 1982 – January 1st 2018, where
satellite observations of SIC from the National Snow and Ice Data Center
(NSIDC; Cavalieri1996) NASA Team algorithm were assimilated into SIS2 every 5
days using the Ensemble Adjustment Kalman Filter (EAKF) approach [Anderson
(2001)], and sea-surface temperatures were nudged towards observations from
version 2 of the Optimum Interpolation Sea-Surface Temperature (OISSTv2) data
set [Reynolds . (2007), Banzon . (2016)] at the model timestep. It should be
noted that SIS2 has a 5-category ice thickness distribution [Bitz . (2001)],
with lower thickness bounds of 0.0, 0.1, 0.3, 0.7, and 1.1 meters. The
(observable) aggregate SIC field is therefore computed in the model as the sum
of the sea ice concentration in each category (SICN), hence SIC
$=\sum_{k=1}^{5}$SICNk. Similarly, we compute the aggregate SIC increment
($\Delta$SIC) as the sum of the analysis increments in each category
($\Delta$SICN). The G23 CNN then uses 5-day mean inputs of state variables and
tendencies corresponding to: SIC, sea-surface temperature (SST), zonal and
meridional components of ice velocities (SIU and SIV, respectively), sea ice
thickness (SIT), net shortwave radiation (SW), ice-surface skin temperature
(TS), sea-surface salinity (SSS), and a land-sea mask, in order to predict
$\Delta$SIC. This prediction of $\Delta$SIC is then passed to a second CNN,
along with SICN, to predict the category concentration increments
$\Delta$SICN. For convenience we refer to these two CNNs as a single network
hereafter.

The implementation of the CNN into SIS2 here is performed in an analogous
manner to DA. Specifically, we run an ensemble forecast of the model for 5
days (e.g., from 00:00 hours UTC on January 1st to 00:00 UTC on January 6th),
where we then generate the corresponding $\Delta$SICN predictions for each
ensemble member, add the predicted $\Delta$SICN fields to the instantaneous
SICN state (i.e., the state at 00:00 UTC on January 6th), and restart the
model for the next 5-day forecast (schematics of this 5-day forecast plus
correction process are shown in Figure 1, although section 2.3 describes this
figure in more detail). It is important to note that we also apply a post-
processing after each correction. For this we follow the Z21 procedure for
updating sea ice variables during DA, which is as follows: first we remove
non-physical values from the updated SICN terms by applying a lower bound of 0
to each category, and then scaling each category by 1/SIC if the updated SIC
is greater than 1. Secondly, in the case where the correction is removing all
sea ice within a given grid cell, we set the corresponding ice and snow
thickness, enthalpy, and ice salinity to 0, and subsequently set the ice-
surface skin temperature to the freezing point of sea water, $-1.8^{\circ}$C.
In the case where the correction is adding sea ice to a given category which
was previously ice-free, we set the thickness of the ice to the mid-point
value within the ice thickness distribution bounds (given as 0.05, 0.2, 0.5,
0.9, 1.3 meters). We then set the salinity, enthalpy and skin temperature of
the ice to 5 psu, $-87576$ J, and $-0.36^{\circ}$C, respectively (conditions
based on an initial liquid fraction of frazil ice of 0.75). We also ensure
that the newly added sea ice contains no overlying snow.

###  2.3 Ice-ocean experiments

We compare four ice-ocean simulations in this study, where each extends for a
5-year period between January 1st 2018 and January 1st 2023\. The initial ice
and ocean conditions for all simulations are based on those from the Z21 DA
experiment, which ended January 1st 2018\. The atmospheric forcing is provided
by JRA55-do reanalysis version 1.5, SSTs are nudged towards OISSTv2
observations using a piston velocity of 4 meters per day, and SSS is nudged to
a seasonal climatology with a piston velocity of 1/6 meters per day. The
experiments are given as follows:

  1. 1.

The free-running model in ice-ocean mode (FREE).

  2. 2.

An extension of the ice-ocean DA experiment (DAZ21), which serves as the
benchmark for this study.

  3. 3.

An ‘out-of-the-box’ implementation of the G23 network (CNNG23), where the
network has been trained offline using all available data from the original DA
experiment. This procedure is highlighted in Figure 1a, where, during
training, the loss function $\mathcal{L}$ minimizes the error between the
network predictions, CNN($\bar{\mathbf{X}}$), and the increment from DA,
$\Delta$SIC(DAZ21). Here $\bar{\mathbf{X}}$ represents the 5-day mean state
variables and tendencies described in section 2.2. After training, this
produces the network CNNG23, which is then implemented over the 2018–2022
period. The reader is referred to G23 for more details of the architecture and
hyperparameters related to the network training process.

  4. 4.

An ‘optimized’ version of the G23 network (CNNopt) where the weights of the
G23 network are refined to improve online performance. For this we use both DA
and CNNs to iteratively augment the training data, and subsequently refine the
network weights after each augmentation iteration. For example, in the first
iteration we run a new ice-ocean simulation between 1982–2017, and in which we
apply a two-step CNN+DA correction every 5 days; first using CNNG23, and then
using DA (see Figure 1b). We then use this 36-year simulation as a new
training data set with which to update the weights of CNNG23, where, during
training, the loss function now minimizes the error between the network
predictions, CNN($\bar{\mathbf{X}}$), and the total model error,
$\Delta$SIC(CNNG23+DA). This procedure is performed for a total of $N=3$
iterations. The network refinement after each augmentation iteration is
performed in an identical way to the offline learning procedure outlined in
G23, except now we only update the weights for 5 epochs after each iteration.

Note that the ‘Online validation’ panels in Figure 1 highlight the simulations
with CNN implementations relative to a simulation which applies the respective
‘perfect’ correction (i.e., the correction which either CNN would produce if
it had 100% prediction accuracy). These simply correspond to the extended DA
experiment for CNNG23 in Figure 1a, and the two-step CNN correction plus DA
for CNNopt in Figure 1b. A comparison of these corrections (increments) is
made in section 3.2.2 in order to establish how the different implementation
configurations manifest within the increments. As a final point to note here,
all results presented in this work are based on ensemble mean fields, and all
simulations are run with a ‘no leap’ calendar, which excludes leap-year days.

Figure 2: Comparison of model biases from ice-ocean simulations over the
2018–2022 period. The first and second rows show the mean SIC biases (model
minus observations), relative to the NSIDC NASA Team observations [DiGirolamo
. (2022)]. The average RMSE of SIC is reported in each panel (RMSE computed
each day over sea ice covered grid cells only, and then averaged over all
days). The bottom two time series plots show pan-Arctic and pan-Antarctic RMSE
of SIC for each simulation.

##  3 Results

###  3.1 Model bias

Figure 2 shows model biases for each of the ice-ocean experiments outlined in
section 2.3. Initially considering the annual-mean spatial bias patterns of
the free-running model, we can see that this simulation is overall positively
biased in both hemispheres, with largest Arctic biases occurring in the east
Atlantic sector (Greenland, Barents, and Kara seas), and largest Antarctic
biases in the Bellingshausen, Amundsen, and Indian Ocean sectors. The daily
SIC root-mean-squared error (RMSE) curves then highlight the seasonal
variation of the model bias, with largest RMSE values in FREE (black curves)
occurring across June–August in the Arctic (22.7%), and December–February in
the Antarctic (28.5%). The average RMSE over the entire simulation period
corresponds to 17.8% and 20.0% in the Arctic and Antarctic, respectively, with
larger errors in summer and smaller errors in winter. As expected, the DA
experiment (DAZ21; red curves) visibly reduces the bias across all seasons,
with average Arctic and Antarctic RMSE reductions relative to FREE of 4.2% and
5.5%, respectively.

Turning to the two CNN correction schemes, the out-of-the-box implementation
(CNNG23; light blue curves) shows systematic improvements relative to FREE,
with average RMSE reductions of 1.9% and 3.9% in the Arctic and Antarctic,
respectively. The modest improvements in the Arctic make it difficult to
identify qualitative differences in the climatology spatial bias plots,
however some improvements can be seen in the east Atlantic sector. This is
also highlighted in the regional Arctic sea ice extent (SIE) time series
(Figure S1), where regions are defined according to Meier2023. On the other
hand, CNNG23 shows visible improvements across much of the Antarctic domain,
with large bias reductions in the Amundsen Sea and Pacific Ocean. The regional
Antarctic SIE time series (Figure S2) also more closely track the DA
experiment throughout the majority of the simulation period, particularly in
the Antarctic growth season. In the melt season however, the simulation shows
a tendency to drift back towards to the free-running model state. Comparing
these results to the optimized CNN implementation (CNNopt; dark blue curves),
we see marked skill improvements. The average RMSE reductions compared to FREE
are 3.9% and 5.9% in the Arctic and Antarctic, respectively, and Figure 2
shows that sizeable improvements have been made in the summer months in both
hemispheres. Furthermore, both pan-Antarctic and regional SIE (Figure S2) are
also considerably improved in the melt season compared to CNNG23, and often
show reduced biases relative to the DA experiment. It is worth noting however
that many of the regional Antarctic SIE time series for CNNopt show visible
imprints of model shock (i.e., large fluctuations in extent, occurring every 5
days). This can occur in DA when there is significant drift between each
assimilation cycle, and the fact that we see this here may suggest that there
is rapid error growth occurring over the space of 5 days in the Antarctic. We
discuss this further in section 3.2.2. In any case, the fact that the CNNopt
experiment, which does not assimilate any observations, has similar errors to
DAZ21 suggests that the DA run is primarily correcting systematic model error
and that CNNopt is successfully capturing these errors.

###  3.2 Understanding online improvements

Between the two CNN models, it is clear that CNNopt is the most desirable
scheme for reducing the free-running model bias. Furthermore, it is also clear
that, relative to CNNG23, the largest gains from CNNopt come in the summer
months. In this section we take a closer look at the performance of each CNN
correction scheme in order to discern how these improvements manifest in both
the spatial error patterns of each simulation, and also in the SIC increments
produced from each scheme.

####  3.2.1 Snapshots

Figure 3 shows example snapshots of summertime model errors in each hemisphere
(see also supplementary movie S1 for snapshots over the 5-year simulation
period). In both hemispheres we can see that FREE (Figures 3a and 3e) contains
large positive errors related to over-estimation of the sea ice edge
(indicated by the positive SIC biases equator-ward of the observed ice edge
contour). The errors pole-ward of the ice edge contour indicate local SIC
errors. While the DA simulation (Figures 3b and 3f) retains some of these
local SIC errors, a significant fraction of the ice edge errors are reduced;
almost halving the RMSE in the Antarctic relative to FREE. CNNG23 (Figures 3c
and 3g) shows some improvements relative to FREE (4.1% and 7.8% RMSE
improvement in the Arctic and Antarctic, respectively), however there are
still considerable ice edge and local SIC errors throughout both the Pacific
sector in the Arctic, and the Atlantic and Pacific sectors in the Antarctic.

For CNNopt (Figures 3d and 3h) there are clear RMSE improvements relative to
both FREE and CNNG23 in both hemispheres, with remarkable improvements in the
Antarctic. It could be argued however that, in the Arctic, the simulated ice
edge position is not much improved in this example. A useful metric to confirm
this is the integrated ice edge error (IIEE; Goessling2016), which computes
the total area for which the ice edge is both over- and under-predicted,
relative to satellite observations. For panels (a–d) in Figure 3, the IIEEs
are given as 1.44, 0.90, 1.32, 1.05 million km2, respectively, which shows
that the sea ice edge from each correction scheme is in better agreement with
the observations than FREE, and that CNNopt does indeed improve over CNNG23 in
this regard. Similarly in the Antarctic panels (e–h), the IIEEs are given as
3.98, 1.38, 3.11, 1.30 million km2, respectively. Here we can see that CNNopt
even shows improved ice edge errors over the DA simulation (see Figure S4 for
IIEE metrics computed over the entire 2018–2022 period). We can also take the
assessment of ice edge errors further by disaggregating the SIC RMSE metric
into grid points that lie pole-ward and equator-ward of the observed ice edge
contour on any given day, in order to assess where the largest improvements
from each correction scheme are manifesting (i.e., whether improvements are
primarily in the ice edge location or SIC within the ice pack). From this
decomposition (Figures S5 and S6) we find that, relative to FREE, in both
hemispheres the largest RMSE reductions from each correction scheme come from
improvements in the ice edge. Furthermore, we find that CNNopt is considerably
reducing the summer ice edge errors relative to CNNG23.

Figure 3: Snapshots of summertime model errors (model minus observations) for
the FREE, DA and CNN ice-ocean simulations. Errors are computed relative to
the NSIDC NASA Team observations [DiGirolamo . (2022)]. RMSE values are
reported in each panel. The black contours mark the observed sea ice edge
position (15% SIC). Figure 4: Comparison of SIC increments produced from
either DA or CNNs during online simulations. The first two rows show spatial
climatologies over 2018–2022. The bottom two panels are then the equivalent
time series, computed as mean fields over Arctic and Antarctic domains.

####  3.2.2 Analysis increments

Figure 4 shows the mean SIC increments for all DA and CNN correction schemes.
The increments here correspond to those which were originally outlined in the
validation panels of Figure 1. Namely, the extended DA experiment (DAZ21), the
out-of-the-box CNN implementation (CNNG23), the optimized version of the G23
network (CNNopt), and each of the corrections from the two-step CNN$+$DA
process, using the optimized network (CNNCNN+DA and DACNN+DA, respectively).

The mean increments from both CNNG23 and CNNopt in Figure 4 show largely
similar spatial patterns in both hemispheres, with CNNopt displaying overall
larger magnitudes. In the Arctic, while both sets of CNN increments show
isolated regions of positive values along the Eurasian coast, they do not
reflect the larger area of mean positive increments seen in DAZ21 across the
East Siberian, Laptev and Kara seas (from the DAZ21 Arctic time series panel
in Figure 4 we can see that these positive increments originate in summer). On
the other hand, the increments from CNNCNN+DA do indeed show mean positive
summer values in these regions. This suggests that the combination of input
variables to the network which are needed to generate these positive
predictions in the Arctic, is only being sampled after the additional DACNN+DA
step. Nonetheless, in section 3.2.1 we have seen that CNNopt yields
significant improvements in Arctic summer SIC errors over CNNG23, which is
primarily coming from improvements in ice edge errors. This is consistent with
the larger magnitude negative increments from CNNopt in regions such as the
Beaufort and Chukchi seas. This may then suggest that the positive summer
increments seen in the DA and CNNCNN+DA corrections are needed to target local
summer SIC errors. Regarding DACNN+DA, we can see that, in both hemispheres,
these corrections are lower in magnitude than DAZ21 on average, which
highlights how the initial correction from CNNCNN+DA is removing a sizeable
component of the model error; leaving less error to correct with DA. This is
particularly the case in the Antarctic, where the daily increments from
DACNN+DA are very close to zero, suggesting that the CNN has effectively
removed the systematic component of the model error in the Antarctic.
Meanwhile in the Arctic, there are still residual systematic summertime errors
associated with under-predicting the positive increments, which DACNN+DA needs
to address.

A natural question then arises as to why DAZ21 is as effective, if not more
effective, at reducing the model bias than CNNopt, even though the increments
from CNNopt are considerably larger in magnitude. For this we turn to a
comparison of the increments from each of the concentration categories
(Figures S6 and S7). For CNNopt, we find that the largest magnitude
corrections are being made to the thinnest ice category, while on the other
hand, DAZ21 makes sizeable corrections to some of the thicker categories. This
therefore means that CNNopt needs to make larger corrections to achieve the
same volume change as DAZ21. Furthermore, the fact that CNNopt is largely
updating the thinnest ice category also explains the model shock seen in the
regional Antarctic SIE time series for CNNopt in Figure S2. In the Pacific
sector in summer for example, the CNN is adding large extents of ice, which
the model is then consistently removing over each 5 day interval. This now
seems conceivable given that the new ice is very thin (5 cm thickness), and
hence would be susceptible to completely melting if advected to grid cells
with sufficiently warm SSTs. Further evidence to support this claim comes from
the fact that this model shock behavior is significantly damped in the
regional sea ice volume time series (Figure S8); which makes sense as the
thinnest ice category will typically contribute less to the regional volume.

##  4 Discussion and conclusions

There is currently much discourse centered around ML-based parameterizations
and/or corrections within climate models, particularly in the context of how
to achieve stable and unbiased simulations after implementation. Many studies
have illustrated how to achieve stability within idealized models, where
parameterizations are learned from high resolution simulations. For example,
by swapping out neural networks for random forests [Yuval O’Gorman (2020),
Watt-Meyer . (2021)], or using online and/or reinforcement learning [Rasp
(2020), Kurz . (2023)].

In this study we have shown that a CNN model which has been trained purely
offline to predict increments from a sea ice DA system (which assimilates real
sea ice observations) can be used ‘out-of-the-box’ to systematically reduce
sea ice biases in a 5-year global ice-ocean simulation, without instabilities
or drift. We have also introduced a data augmentation approach to optimize the
offline-trained CNN, which significantly improves online generalization in
both hemispheres; particularly in terms of reducing sea ice edge errors in
summer. This augmentation approach is performed by iteratively generating new
simulations in which corrections are applied from both the current iteration
of the ML model, as well as DA. Each iteration of the augmentation procedure
therefore provides a new training data set with which to refine the CNN
weights from the previous iteration. While, in theory, this procedure could be
repeated to convergence, we opted for $N=3$ iterations in this study due to
computational expense. It is likely however that continued iterations would
yield further improvements, particularly in the Arctic summer, given that
sizeable gains were made between each of the three iterations here (see Figure
S9). We hypothesise that the improvements from this augmentation procedure are
a result of exposing the network to input variables which contain information
about how the model trajectory evolves after implementing the CNN (as opposed
to training purely offline where the inputs have no feedback with the CNN).

Interestingly, we find that, relative to the original DA experiment, the
climatological sea ice biases associated with the simulation which uses this
‘optimized’ network are actually modestly improved in the Antarctic (Figure
2). This is understandable when we consider that the target variable during
each iteration of the network refinement step is no longer the increment from
original DA experiment, but rather the sum of the increments from the two-step
CNN$+$DA experiment (recall Figure 1b). Therefore, the model bias from the
original DA experiment should not be seen as the lower limit on what is
achievable with the CNN. We can see this in Figure S10, where the bias of the
simulation which applies this two-step CNN$+$DA procedure is indeed
systematically lower than both the original DA experiment and the optimized
CNN. This leaves exciting avenues for future work relating to improved initial
conditions for numerical prediction. For seasonal predictions with the GFDL
SPEAR model for example, initial conditions for the ice and ocean [Y. Zhang .
(2022), Lu . (2020)] are based on DA via Ensemble Kalman filters. The Ensemble
Kalman filter is not formally designed to correct for systematic model error,
and so this two-step CNN$+$DA procedure could be a way to generate more
accurate initial conditions (see Figure S10) with the CNN and DA fixing the
systematic and random components of the errors, respectively. Indeed while
this is similar to a weak constraint 4-D variational DA approach [Wergen
(1992), Zupanski (1993), Trémolet (2007)], the CNN has computational
advantages (once trained) in that it does not require the construction of an
adjoint model [Bonavita Laloyaux (2020), Laloyaux . (2022)].

Further avenues for future work also include the use of the optimized CNN for
making bias corrections to real-time seasonal sea ice forecasts, or sea ice
projections on climate timescales. For seasonal prediction, the methodology
would follow that which has been presented in this study, except the CNN would
be applied within the fully coupled SPEAR model. This would however require
several considerations. For example, addressing the issue of model shock seen
in the Antarctic (which could potentially be reduced by increasing the
frequency of the CNN corrections), and also assessing generalization of the
CNN to the fully coupled model, which includes new interactive feedbacks with
an atmospheric model. Looking also to longer term climate projections, G23
discussed this in the context of implementing the CNN as a sea ice model
parameterization. This would then require further considerations of how to
appropriately conserve mass, heat and salt when adding/removing sea ice from
the ocean.

##  5 Open Research

All data for training each CNN are openly available [Gregory (2023)], along
with auxiliary data such as the optimized CNN weights and standardization
statistics. Python code to pre-process the input data and train the CNNs is
also available at the same location.

###### Acknowledgements.

William Gregory, Mitchell Bushuk, Alistair Adcroft and Laure Zanna received
M2LInES research funding by the generosity of Eric and Wendy Schmidt by
recommendation of the Schmidt Futures program. This work was also
intellectually supported by various other members of the M2LInES project, as
well as being supported through the provisions of computational resources from
the National Oceanic and Atmospheric Administration (NOAA) Geophysical Fluid
Dynamics Laboratory (GFDL). We also thank Theresa Morrison and Feiyu Lu for
their invaluable feedback on this article.

## References

  * Adcroft . (2019) Adcroft2019Adcroft, A., Anderson, W., Balaji, V., Blanton, C., Bushuk, M., Dufour, CO.others  2019\.  The GFDL global ocean and sea ice model OM4.0: Model description and simulation features The GFDL global ocean and sea ice model OM4.0: Model description and simulation features.  Journal of Advances in Modeling Earth Systems11103167–3211.  https://doi.org/10.1029/2019MS001726
  * Anderson (2001) Anderson2001Anderson, J.  2001\.  An ensemble adjustment Kalman filter for data assimilation An ensemble adjustment Kalman filter for data assimilation.  Monthly weather review129122884–2903.  https://doi.org/10.1175/1520-0493(2001)129¡2884:AEAKFF¿2.0.CO;2
  * Banzon . (2016) Banzon2016Banzon, V., Smith, TM., Chin, TM., Liu, C. Hankins, W.  2016\.  A long-term record of blended satellite and in situ sea-surface temperature for climate monitoring, modeling and environmental studies A long-term record of blended satellite and in situ sea-surface temperature for climate monitoring, modeling and environmental studies.  Earth System Science Data81165–176.  https://doi.org/10.5194/essd-8-165-2016
  * Bitz . (2001) Bitz2001Bitz, CM., Holland, MM., Weaver, AJ. Eby, M.  2001\.  Simulating the ice-thickness distribution in a coupled climate model Simulating the ice-thickness distribution in a coupled climate model.  Journal of Geophysical Research: Oceans106C22441–2463.  https://doi.org/10.1029/1999JC000113
  * Bonavita Laloyaux (2020) Bonavita2020Bonavita, M. Laloyaux, P.  2020\.  Machine learning for model error inference and correction Machine learning for model error inference and correction.  Journal of Advances in Modeling Earth Systems1212e2020MS002232.  https://doi.org/10.1029/2020MS002232
  * Brajard . (2021) Brajard2021Brajard, J., Carrassi, A., Bocquet, M. Bertino, L.  2021\.  Combining data assimilation and machine learning to infer unresolved scale parametrization Combining data assimilation and machine learning to infer unresolved scale parametrization.  Philosophical Transactions of the Royal Society A379219420200086.  https://doi.org/10.1098/rsta.2020.0086
  * Brenowitz . (2020) Brenowitz2020Brenowitz, ND., Beucler, T., Pritchard, M. Bretherton, CS.  2020\.  Interpreting and stabilizing machine-learning parametrizations of convection Interpreting and stabilizing machine-learning parametrizations of convection.  Journal of the Atmospheric Sciences77124357–4375.  https://doi.org/10.1175/JAS-D-20-0082.1
  * Brenowitz Bretherton (2018) Brenowitz2018Brenowitz, ND. Bretherton, CS.  2018\.  Prognostic validation of a neural network unified physics parameterization Prognostic validation of a neural network unified physics parameterization.  Geophysical Research Letters45126289–6298.  https://doi.org/10.1029/2018GL078510
  * Cavalieri . (1996) Cavalieri1996Cavalieri, DJ., Parkinson, CL., Gloersen, P. Zwally, HJ.  1996\.  Sea ice concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS passive microwave data, version 1. NASA Natl. Snow and Ice Data Cent Sea ice concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS passive microwave data, version 1. NASA Natl. Snow and Ice Data Cent.  Distrib. Active Arch. Cent., Boulder, Colo..  https://doi.org/10.5067/8GQ8LZQVL0VL
  * Chen . (2022) Chen2022Chen, TC., Penny, SG., Whitaker, JS., Frolov, S., Pincus, R. Tulich, S.  2022\.  Correcting systematic and state-dependent errors in the NOAA FV3-GFS using neural networks Correcting systematic and state-dependent errors in the NOAA FV3-GFS using neural networks.  Journal of Advances in Modeling Earth Systems1411.  https://doi.org/10.1029/2022MS003309
  * Delworth . (2020) Delworth2020Delworth, TL., Cooke, WF., Adcroft, A., Bushuk, M., Chen, JH., Dunne, KA.others  2020\.  SPEAR: The next generation GFDL modeling system for seasonal to multidecadal prediction and projection SPEAR: The next generation GFDL modeling system for seasonal to multidecadal prediction and projection.  Journal of Advances in Modeling Earth Systems123e2019MS001895.  https://doi.org/10.1029/2019MS001895
  * DiGirolamo . (2022) DiGirolamo2022DiGirolamo, N., Parkinson, C., Cavalieri, D., Gloersen, P. Zwally, H.  2022\.  Sea ice concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS passive microwave data, version 2. Sea ice concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS passive microwave data, version 2.  National Snow and Ice Data Center.  https://doi.org/10.5067/MPYG15WAA4WX
  * Farchi . (2021) Farchi2021Farchi, A., Laloyaux, P., Bonavita, M. Bocquet, M.  2021\.  Using machine learning to correct model error in data assimilation and forecast applications Using machine learning to correct model error in data assimilation and forecast applications.  Quarterly Journal of the Royal Meteorological Society1477393067–3084.  https://doi.org/10.1002/qj.4116
  * Finn . (2023) Finn2023Finn, T., Durand, C., Farchi, A., Bocquet, M., Chen, Y., Carrassi, A. Dansereau, V.  2023\.  Deep learning of subgrid-scale parametrisations for short-term forecasting of sea-ice dynamics with a Maxwell-Elasto-Brittle rheology Deep learning of subgrid-scale parametrisations for short-term forecasting of sea-ice dynamics with a maxwell-elasto-brittle rheology.  EGUsphere2022-1342.  https://doi.org/10.5194/egusphere-2022-1342
  * Frezat . (2022) Frezat2022Frezat, H., Le Sommer, J., Fablet, R., Balarac, G. Lguensat, R.  2022\.  A posteriori learning for quasi-geostrophic turbulence parametrization A posteriori learning for quasi-geostrophic turbulence parametrization.  Journal of Advances in Modeling Earth Systems1411e2022MS003124.  https://doi.org/10.1029/2022MS003124
  * Gentine . (2018) Gentine2018Gentine, P., Pritchard, M., Rasp, S., Reinaudi, G. Yacalis, G.  2018\.  Could machine learning break the convection parameterization deadlock? Could machine learning break the convection parameterization deadlock?  Geophysical Research Letters45115742–5751.  https://doi.org/10.1029/2018GL078202
  * Goessling . (2016) Goessling2016Goessling, HF., Tietsche, S., Day, JJ., Hawkins, E. Jung, T.  2016\.  Predictability of the Arctic sea ice edge Predictability of the arctic sea ice edge.  Geophysical Research Letters4341642–1650.  https://doi.org/10.1002/2015GL067232
  * Gregory (2023) Gregory2023bGregory, W.  2023\.  Sea ice DA-ML: April 11, 2023 release (version 1.0.0) [Dataset, Software]. Sea ice DA-ML: April 11, 2023 release (version 1.0.0) [Dataset, Software].  https://doi.org/10.5281/zenodo.7818178
  * Gregory . (2023) Gregory2023Gregory, W., Bushuk, M., Adcroft, A., Zhang, Y. Zanna, L.  2023\.  Deep learning of systematic sea ice model errors from data assimilation increments Deep learning of systematic sea ice model errors from data assimilation increments.  Journal of Advances in Modeling Earth Systems15e2023MS003757.  https://doi.org/10.1029/2023MS003757
  * He . (2023) He2023He, Z., Brajard, J., Wang, Y., Wang, X. Shen, Z.  2023\.  Improve dynamical climate prediction with machine learning Improve dynamical climate prediction with machine learning.  ESS Open Archive.  https://doi.org/10.22541/essoar.167898499.91486328/v1
  * Held . (2019) Held2019Held, I., Guo, H., Adcroft, A., Dunne, J., Horowitz, L., Krasting, J.others  2019\.  Structure and performance of GFDL’s CM4. 0 climate model Structure and performance of GFDL’s CM4. 0 climate model.  Journal of Advances in Modeling Earth Systems11113691–3727.  https://doi.org/10.1029/2019MS001829
  * Kurz . (2023) Kurz2023Kurz, M., Offenhäuser, P. Beck, A.  2023\.  Deep reinforcement learning for turbulence modeling in large eddy simulations Deep reinforcement learning for turbulence modeling in large eddy simulations.  International Journal of Heat and Fluid Flow99109094.  https://doi.org/10.1016/j.ijheatfluidflow.2022.109094
  * Laloyaux . (2022) Laloyaux2022Laloyaux, P., Kurth, T., Dueben, PD. Hall, D.  2022\.  Deep learning to estimate model biases in an operational NWP assimilation system Deep learning to estimate model biases in an operational NWP assimilation system.  Journal of Advances in Modeling Earth Systemse2022MS003016.  https://doi.org/10.1029/2022MS003016
  * Lu . (2020) Lu2020Lu, F., Harrison, MJ., Rosati, A., Delworth, TL., Yang, X., Cooke, WF.others  2020\.  GFDL’s SPEAR seasonal prediction system: Initialization and ocean tendency adjustment (OTA) for coupled model predictions GFDL’s SPEAR seasonal prediction system: Initialization and ocean tendency adjustment (OTA) for coupled model predictions.  Journal of Advances in Modeling Earth Systems1212e2020MS002149.  https://doi.org/10.1029/2020MS002149
  * Meier Stewart (2023) Meier2023Meier, WN. Stewart, JS.  2023\.  Arctic and Antarctic regional masks for sea ice and related data products, version 1. Arctic and Antarctic regional masks for sea ice and related data products, version 1.  National Snow and Ice Data Center.  https://doi.org/10.5067/CYW3O8ZUNIWC
  * Mojgani . (2022) Mojgani2022Mojgani, R., Chattopadhyay, A. Hassanzadeh, P.  2022\.  Discovery of interpretable structural model errors by combining Bayesian sparse regression and data assimilation: A chaotic Kuramoto–Sivashinsky test case Discovery of interpretable structural model errors by combining Bayesian sparse regression and data assimilation: A chaotic Kuramoto–Sivashinsky test case.  Chaos: An Interdisciplinary Journal of Nonlinear Science326061105.  https://doi.org/10.1063/5.0091282
  * O’Gorman Dwyer (2018) OGorman2018O’Gorman, PA. Dwyer, JG.  2018\.  Using machine learning to parameterize moist convection: Potential for modeling of climate, climate change, and extreme events Using machine learning to parameterize moist convection: Potential for modeling of climate, climate change, and extreme events.  Journal of Advances in Modeling Earth Systems10102548–2563.  https://doi.org/10.1029/2018MS001351
  * Ott . (2020) Ott2020Ott, J., Pritchard, M., Best, N., Linstead, E., Curcic, M. Baldi, P.  2020\.  A Fortran-Keras deep learning bridge for scientific computing A Fortran-Keras deep learning bridge for scientific computing.  Scientific Programming2020.  https://doi.org/10.1155/2020/8888811
  * Rasp (2020) Rasp2020Rasp, S.  2020\.  Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: general algorithms and Lorenz 96 case study (v1. 0) Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: general algorithms and Lorenz 96 case study (v1. 0).  Geoscientific Model Development1352185–2196.  https://doi.org/10.5194/gmd-13-2185-2020
  * Rasp . (2018) Rasp2018Rasp, S., Pritchard, MS. Gentine, P.  2018\.  Deep learning to represent subgrid processes in climate models Deep learning to represent subgrid processes in climate models.  Proceedings of the National Academy of Sciences115399684–9689.  https://doi.org/10.1073/pnas.1810286115
  * Reynolds . (2007) Reynolds2007Reynolds, RW., Smith, TM., Liu, C., Chelton, DB., Casey, KS. Schlax, MG.  2007\.  Daily high-resolution-blended analyses for sea surface temperature Daily high-resolution-blended analyses for sea surface temperature.  Journal of climate20225473–5496.  https://doi.org/10.1175/2007JCLI1824.1
  * Ross . (2023) Ross2023Ross, A., Li, Z., Perezhogin, P., Fernandez-Granda, C. Zanna, L.  2023\.  Benchmarking of machine learning ocean subgrid parameterizations in an idealized model Benchmarking of machine learning ocean subgrid parameterizations in an idealized model.  Journal of Advances in Modeling Earth Systems151e2022MS003258.  https://doi.org/10.1029/2022MS003258
  * Sane . (2023) Sane2023Sane, A., Reichl, BG., Adcroft, A. Zanna, L.  2023\.  Parameterizing Vertical Mixing Coefficients in the Ocean Surface Boundary Layer using Neural Networks Parameterizing vertical mixing coefficients in the ocean surface boundary layer using neural networks.  Journal of Advances in Modeling Earth Systems15e2023MS003890.  https://doi.org/10.1029/2023MS003890
  * Trémolet (2007) Tremolet2007Trémolet, Y.  2007\.  Incremental 4D-Var convergence study Incremental 4d-Var convergence study.  Tellus A: Dynamic Meteorology and Oceanography595706–718.  https://doi.org/10.1111/j.1600-0870.2007.00271.x
  * Tsujino . (2018) Tsujino2018Tsujino, H., Urakawa, S., Nakano, H., Small, RJ., Kim, WM., Yeager, SG.others  2018\.  JRA-55 based surface dataset for driving ocean–sea-ice models (JRA55-do) JRA-55 based surface dataset for driving ocean–sea-ice models (JRA55-do).  Ocean Modelling13079–139.  https://doi.org/10.1016/j.ocemod.2018.07.002
  * Watt-Meyer . (2021) Watt2021Watt-Meyer, O., Brenowitz, ND., Clark, SK., Henn, B., Kwa, A., McGibbon, J.Bretherton, CS.  2021\.  Correcting weather and climate models by machine learning nudged historical simulations Correcting weather and climate models by machine learning nudged historical simulations.  Geophysical Research Letters4815e2021GL092555.  https://doi.org/10.1029/2021GL092555
  * Wergen (1992) Wergen1992Wergen, W.  1992\.  The effect of model errors in variational assimilation The effect of model errors in variational assimilation.  Tellus A444297–313.  https://doi.org/10.1034/j.1600-0870.1992.t01-3-00002.x
  * Yuval O’Gorman (2020) Yuval2020Yuval, J. O’Gorman, PA.  2020\.  Stable machine-learning parameterization of subgrid processes for climate modeling at a range of resolutions Stable machine-learning parameterization of subgrid processes for climate modeling at a range of resolutions.  Nature communications1111–10.  https://doi.org/10.1038/s41467-020-17142-3
  * Zanna Bolton (2020) Zanna2020Zanna, L. Bolton, T.  2020\.  Data-driven equation discovery of ocean mesoscale closures Data-driven equation discovery of ocean mesoscale closures.  Geophysical Research Letters4717e2020GL088376.  https://doi.org/10.1029/2020GL088376
  * C. Zhang . (2023) Zhang2023Zhang, C., Perezhogin, P., Gultekin, C., Adcroft, A., Fernandez-Granda, C. L, Z.  2023\.  Implementation and evaluation of a machine learned mesoscale eddy parameterization into a numerical ocean circulation model Implementation and evaluation of a machine learned mesoscale eddy parameterization into a numerical ocean circulation model.  Journal of Advances in Modeling Earth Systems15e2023MS003697.  https://doi.org/10.1029/2023MS003697
  * Y. Zhang . (2022) Zhang2022Zhang, Y., Bushuk, M., Winton, M., Hurlin, B., Delworth, T., Harrison, M.Yang, X.  2022\.  Subseasonal-to-seasonal Arctic sea ice forecast skill improvement from sea ice concentration assimilation Subseasonal-to-seasonal Arctic sea ice forecast skill improvement from sea ice concentration assimilation.  Journal of Climate1–48.  https://doi.org/10.1175/JCLI-D-21-0548.1
  * Y. Zhang . (2021) Zhang2021Zhang, Y., Bushuk, M., Winton, M., Hurlin, B., Yang, X., Delworth, T. Jia, L.  2021\.  Assimilation of satellite-retrieved sea ice concentration and prospects for september predictions of Arctic sea ice Assimilation of satellite-retrieved sea ice concentration and prospects for september predictions of Arctic sea ice.  Journal of Climate3462107–2126.  https://doi.org/10.1175/JCLI-D-20-0469.1
  * Zupanski (1993) Zupanski1993Zupanski, M.  1993\.  Regional four-dimensional variational data assimilation in a quasi-operational forecasting environment Regional four-dimensional variational data assimilation in a quasi-operational forecasting environment.  Monthly weather review12182396–2408.  https://doi.org/10.1175/1520-0493(1993)121¡2396:RFDVDA¿2.0.CO;2

Atmospheric and Oceanic Sciences Program, Princeton University, NJ, USA
Geophysical Fluid Dynamics Laboratory, NOAA, Princeton, NJ, USA Courant
Institute of Mathematical Sciences, New York University, New York, NY, USA

Figure S1: Pan-Arctic and regional daily SIE for ice-ocean experiments over
the 2018–2022 period. The experiments correspond to those described in section
2.3 of the main article. The observations correspond to SIE computed from the
NSIDC NASA Team SIC data set [DiGirolamo . (2022)] The sea ice mask used to
compute regional extents was also derived from the NSIDC product [Meier
Stewart (2023)]. Figure S2: Pan-Antarctic and regional daily SIE for ice-ocean
experiments over the 2018–2022 period. The experiments correspond to those
described in section 2.3 of the main article. The observations correspond to
SIE computed from the NSIDC NASA Team SIC data set [DiGirolamo . (2022)] The
sea ice mask used to compute regional extents was also derived from the NSIDC
product [Meier Stewart (2023)]. Figure S3: Integrated ice edge error (IIEE)
for the various ice-ocean simulations over the 2018–2022 period. IIEE
[Goessling . (2016)] quantifies error in the sea ice edge position as the sum
of grid cell areas which are either over-predicted or under-predicted. IIEE is
computed relative to NSIDC NASA Team SIC observations [DiGirolamo . (2022)].
Figure S4: Disaggregation of Arctic RMSE of SIC for CNN simulations, into
‘local ice errors’ (defined as grid points where the observed SIC is greater
than or equal to 15%), and ‘ice edge errors’ (defined as grid points where the
observed SIC is less than 15%). Observations correspond to the NSIDC NASA Team
data set [DiGirolamo . (2022)]. Figure S5: Disaggregation of Antarctic RMSE of
SIC for CNN simulations, into ‘local ice errors’ (defined as grid points where
the observed SIC is greater than or equal to 15%), and ‘ice edge errors’
(defined as grid points where the observed SIC is less than 15%). Observations
correspond to the NSIDC NASA Team data set [DiGirolamo . (2022)]. Figure S6:
Comparison of category SIC increments in the Arctic, produced from either DA
or CNNs during online simulations. Shown as climatologies over 2018–2022.
Lower thickness bounds of categories 1 through 5 are: 0.0, 0.1, 0.3, 0.7, and
1.1 meters, respectively. Figure S7: Comparison of category SIC increments in
the Antarctic, produced from either DA or CNNs during online simulations.
Shown as climatologies over 2018–2022. Lower thickness bounds of categories 1
through 5 are: 0.0, 0.1, 0.3, 0.7, and 1.1 meters, respectively. Figure S8:
Pan-Antarctic and regional sea ice volume for ice-ocean simulations over the
2018–2022 period. The experiments correspond to those described in section 2.3
of the main article. The sea ice mask used to compute regional volume was
derived from the NSIDC product [Meier Stewart (2023)]. Figure S9: RMSE of SIC
for FREE and CNN simulations over 2018–2022. Here the RMSE is shown for each
iteration (IT) of the optimized CNN. Figure S10: RMSE of SIC for FREE, DA, and
CNN+DA simulations over 2018–2022.


# Efficient Reservoir Computing using Field Programmable Gate Array and
Electro-optic Modulation

Prajnesh Kumar  1 Mingwei Jin  1  Ting Bu1  Santosh Kumar1 and Yu-Ping
Huang1,2  1Physics Department and the Center for Quantum Science and
Engineering, Stevens Institute of Technology, Hoboken, NJ, 07030, USA.
2Department of Electrical and Computer Engineering, Stevens Institute of
Technology, Hoboken, NJ, 07030, USA<EMAIL_ADDRESS>
###### Abstract

We experimentally demonstrate a hybrid reservoir computing system consisting
of an electro-optic modulator and field programmable gate array (FPGA). It
implements delay lines and filters digitally for flexible dynamics and high
connectivity, while supporting a large number of reservoir nodes. To evaluate
the system’s performance and versatility, three benchmark tests are performed.
The first is the $\bm{10^{th}}$ order Nonlinear Auto-Regressive Moving Average
test (NARMA-10), where the predictions of 1000 and 25,000 steps yield
impressively low normalized root mean square errors (NRMSE’s) of 0.142 and
0.148, respectively. Such accurate predictions over into the far future speak
to its capability of large sample size processing, as enabled by the present
hybrid design. The second is the Santa Fe laser data prediction, where a
normalized mean square error (NMSE) of 6.73$\bm{\times 10^{-3}}$ is
demonstrated. The third is the isolate spoken digit recognition, with a word
error rate close to 0.34%. Accurate, versatile, flexibly reconfigurable, and
capable of long-term prediction, this reservoir computing system could find a
wealth of impactful applications in real-time information processing, weather
forecasting, and financial analysis.

††journal: osac††articletype: Research Article

##  1 Introduction

Modern computers based on von Neumann architectures have been designed for
digital information processing and generic computational tasks there upon.
With more and more problems being identified as computationally complex and/or
excessively time consuming, alternative computational paradigms are actively
explored to solve these problems as they increasingly emerge [1, 2, 3, 4, 5].
To that end, a promising approach is to develop brain-inspired architectures
for information processing, including a variety of neural networks [6, 7, 8,
9]. Yet, those architectures are still implemented through computer
simulations on Von Neumann computers, thus fundamentally subject to the
latter’s limitations in speed, parallelism, etc.

Recently, there has been increasing interest on artificial neural networks
using optics, leveraging its remarkable speed, multiplexing capability, and
little heat deposition [10, 11]. Usually, these optical neural networks are
wholly trained with a known data set to optimize their connectivity and
parameters through nonlinear layers [12, 13]. However, such training is
usually energy and time consuming, and its efficiency varies by the complexity
of task, the size of the network, the nonlinearity and connectivity between
the nodes. Also, it is conducted offline using digital computers, thus failing
to account for the uncertainties and fluctuations that are inevitable in any
optics systems, especially when the optical networks are complex and involve
many free-space parts [14, 15].

In an effort to address the above difficulties, the idea of reservoir
computing (RC) was proposed and widely explored [16, 17, 18, 19, 20]. RC
origins from liquid-state machine (LSM) [21] and echo state networks (ESN)
[22]. Its realizations are generally composed of three parts: an input layer,
a reservoir layer, and an output layer. The input signals are first fed into
the input layer, then mapped to the reservoir layer, which contains a large
number of interconnected nonlinear nodes and performs nonlinear
transformation. Afterwards, response is readout with linear weighted sum of
the reservoir states in the output layer. Unlike other kinds of recurrent
neural networks that are notoriously hard to be trained, here only the output
weights are needed to be optimized, making the training much simpler. It is
this advantage and its success on lots of time dependent tasks, such as
chaotic time series prediction and speech recognition [23, 24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35], that draw much attention across different
application areas.

The first RC system was implemented using analog electronics with a single
nonlinear node and delayed feedback [17]. Since then, digital electronics are
increasingly adopted [36, 37, 38]. Among them, field programmable gate array
(FPGA) electronics can make the system compact, stable, low-cost and easily
configurable with various commercial systems for real-time information
processing. More recently, hybrid opto-electronic feedback systems with FPGA’s
have reformulated the development of a coherent Ising machine for solving
computationally hard optimization problems [39] and to build a large networks
of identical nodes with arbitrary topology for cluster synchronization,
chimera states [40] and laminar chaos [41].

Here, we add to the exciting progress made in the RC field and experimentally
demonstrate a fully-packaged opto-electronic RC system consisting of a Mach-
Zehnder electro-optic modulator (EOM) and a FPGA circuit. In our design, both
the delay line and filters are implemented digitally within FPGA, which
renders the whole system compact and immune to optical drifts and noise,
especially compared to fiber-optical realizations [42, 43, 44]. Leveraging the
filters inside the FPGA, we are able to achieve more dynamics and connections
between reservoir nodes. To characterize the system, we run three benchmark
tasks: the $10^{th}$ order Non-Linear Auto-Regressive Moving Average test
(NARMA-10), the Santa Fe laser data prediction, and the isolate spoken digit
recognition. All exhibit exceptionally high performance, which indicates the
robustness and versatility of this RC system.

The remainder of the paper is organized as follows. Section 2 discusses the
basics theory of opto-electronic RC system using delay feedback. In Section 3,
we show the experimental setup which covers detailed block level FPGA
implementations. Subsequently, Section 4 shows the experimental results to
discuss three benchmark tests for evaluating the performance of RC system and
finally, Section 5 concludes this paper.

##  2 Theory

Figure 1: (a) Conventional reservoir computing (RC) model and (b) time-delay
based RC model.

Figure 1(a) shows the conventional RC model which consists of an input, a
reservoir, and an output layer. The input layer consists of input vector
u($t$) of length $K$ which are fed into the reservoir via fixed but random
input weighted connection. These weights will scale the input differently for
different reservoir nodes. The reservoir layer contains a large number of
recurrent and randomly interconnected nonlinear nodes. Reservoir layer non-
linearly projects the inputs to high-dimensional state space. The dynamics of
the reservoir also exhibit a fading memory where the current reservoir state
is influenced only by the recent past. The dynamics of states in the reservoir
layer is given by

$x_{i}(t)=f_{NL}\Big{(}\sum_{k}w_{ik}x_{k}(t-\tau)+\sum_{j}M_{ij}u_{j}(t)\Big{)},$
(1)

where $f_{NL}$ is the non-linear function, $x_{i}(t)$ is the reservoir state
of $i^{th}$ node at time $t$, $w_{ik}$ is the nodes inter-connection matrix,
$M_{ij}$ is the input weight matrix and $u_{j}(t)$ is the $j^{th}$ input.

At the output layer, the response is readout by linear weighted sum of node
states which described as

$\hat{y}_{j}(t)=\sum_{i}W_{ij}x_{i}(t),$ (2)

where $W_{ij}$ being the output weights. The output weights $W_{ij}$ are
updated after training to make the outputs $\hat{y}_{j}(t)$ as close as
possible to the target values $y_{j}(t)$ using optimization methods such as
linear regression, ridge regression etc.

Above RC architecture requires very large number of interconnected physical
nodes. In contrast, the time-delay based RC uses virtual nodes that are
temporally spaced and has only one nonlinear node as shown in Fig. 1(b). The
non-linearity is implemented electro-optically, using an EOM. The input is
preprocessed before injecting into nonlinear node, this procedure is referred
as masking. The preprocessed input signal is time-multiplexed and injected
serially into reservoir. The input vector after preprocessing will resides in
total delay time of $\tau$. The temporal spacing between N nodes is given by
$\theta=\tau/N$ as shown in Fig. 1(b).

Figure 2: Schematic representation of digital opto-electronic RC

Figure 2 shows the schematic of our digital opto-electronic RC implementation
of delay feedback system. The optical power at the output of EOM is given by a
transfer function

$P(t)\sim\sin^{2}[\pi v(t)/(2V_{\pi})+\phi],$ (3)

where $v(t)$ is the RF signal applied to EOM, $V_{\pi}$ is the $\pi$\- shift
voltage and $\phi$ is the bias offset.

The photo-detector generates voltage in response to the applied laser power
which is given by

$V_{det}(t)=\alpha\eta R_{t}P(t),$ (4)

where $\alpha$ is the total insertion loss of the EOM, $\eta$ is the
responsivity of the detector and $R_{t}$ is transimpedance gain of the
amplifier. The RF input to EOM is

$v(t)=G[s(t)+\gamma u(t)],$ (5)

where $G$ is the overall forward gain of the system,$\gamma$ is the input
scaling factor and $u(t)$ is the input to the reservoir. $s(t)$ is the output
of the filter after delay, which is convolution of the delayed detector output
and filter impulse response function $h(t)$, as

$\begin{split}s(t)&=h(t)*V_{det}(t-\tau),\\\ \end{split}$ (6)

with $\tau$ the total delay.

Taking all into account, the equation of the states reads

$\begin{split}x(t)&=\beta
h(t)*\sin^{2}\left[x(t-\tau)+\gamma^{\prime}u(t-\tau)+\phi\right]\\\
\end{split}$ (7)

with $x(t)=\pi Gs(t)/2V_{\pi}$ , $\beta=\pi G\eta
R_{t}\alpha|E_{0}|^{2}/2V_{\pi}$, and $\gamma^{\prime}=G\gamma/2V_{\pi}$. It
is then discretized as

$\begin{split}x[k]&=\beta
h[k]*\sin^{2}\left(x[k-N]+\gamma^{\prime}u[k-N]+\phi\right),\\\ \end{split}$
(8)

where $k=t/dt$ is discrete sample obtained with sampling period $dt$ and
$\tau$ is integer multiple $N$ of $dt$. For digital filter, $h[k]$ has a
finite length of $M$ with $M\leq N$. Thus, Equation (8) is in an explicit form
of

$\begin{split}x[k]=\beta\sum_{j=0}^{M-1}h[j]\sin^{2}\left(x[k-N-j]+\gamma^{\prime}u[k-N-j]+\phi\right).\end{split}$
(9)

From Eq. (9), we can see that the states are coupled using filter
coefficients. Digital filter gives the flexibility to implement different
network topology [42].

##  3 Experimental Setup

The block diagram of our experiment setup is depicted in Figure 3(a). A fiber-
coupled laser diode (JDS Uniphase CQF975/58) provides a continuous wave beam
at $\sim$1550 nm, which is passed through a fiber polarization controller
(FPC) before coupling to an EOM (Lucent Technologies, 2623NA). The EOM takes
RF input signal from the arbitrary waveform generated using RF Digital to
Analog converter (DAC) and modulates the laser beam intensity. The modulated
intensity is converted to electrical signal using a photodetector (Thorlabs,
DET10C2), whose output is digitized using an analog to digital converter (ADC)
with 1 Msps sampling rate controlled by FPGA. A Personal computer (PC)
provides the preprocessed data via Ethernet interface for this RC system.

Figure 3:  (a) Experimental setup for the present opto-electronic RC and (b)
the block diagram of the hardware implementation using FPGA board, and
ADC/DAC.

Figure 3(b) shows the FPGA blocks implemented on a Zedboard (Zynq-7000)
development platform. The RC logic is implemented using Verilog programming
language. The FPGA interfaces to ADC (Maxim Integrated MAX11198, 16bits) and
RF DAC (Analog Devices D5541A,16bits) using Serial Peripheral Interface. The
PC preprocess the input signal using mask and input scaling factor, which is
then sent to a Zedboard via Ethernet interface. Processing System (PS) in
Zedboard shares data with Direct Memory Access (DMA) to stream the input to
the FPGA logic. The feedback signal from filter block and the streamed input
are added and passed through programmable gain and offset block before
converting to electrical signal using DAC. Here, the gain and offset are key
parameters for tuning RC system performance. Similarly streaming data from ADC
is fed into programmable delay block with delay limited to 1000 units, each
delay representing 1 virtual node spacing. The delayed signal is then passed
to finite digital bandpass filter which has 400 filter taps. The filtered data
implements the Eq. (9) which has state information. PS will then transmit this
data to PC via Ethernet interface. Calibration DAC (Maxim Integrated MAX5316,
16bits) is used to compensate for the bias drift of EOM.

##  4 Results

###  4.1 NARMA10

NARMA is one of the most widely used benchmarks in RC[45]. NARMA10 is a
discrete time nonlinear task with $10^{th}$ order lag. Series $y_{k+1}$ is
generated through a recursive formula

$y_{k+1}=0.3y_{k}+0.05y_{k}\sum_{i=0}^{9}y_{k-i}+1.5u_{k}u_{k-9}+0.1.$ (10)

The input $u_{k}$ is drawn from a uniform distribution in the interval
$[0,0.5]$. Due to non-linearity and long time lag, NARMA10 poses a challenge
for any computation system.

To characterize the performance of the RC, a normalized root mean squared
error (NRMSE) between the target and predicted value is calculated as

$\textrm{NRMSE}=\sqrt{\frac{1}{m}\frac{\sum_{k=0}^{m}(\hat{y}_{k}-y_{k})^{2}}{\sigma^{2}(y_{k})}},$
(11)

where $y_{k}$ is the target, $\hat{y}_{k}$ is the prediction, $m$ is the total
number of samples in the target and $\sigma$ denotes the standard deviation of
the target.

By sweeping the system parameters, the optimum operating point is identified
at gain $G=0.58$, input scaling $\gamma=0.5$, bias $\phi=0.1\pi$, total of
$N=400$ virtual nodes. The results of benchmark are tabulated in Table 1.
Figure 4 plots the target vs prediction for 2000 samples. Also plotted is the
residue $\mathrm{\bf R}$, defined as the difference between target and
estimate normalized by the mean of target

$\mathrm{\bf R}[k]=m(y_{k}-\hat{y}_{k})/\sum_{k}y(k).$ (12)

As seen, the present RC performs remarkably well over long term prediction.
For 25000 samples with training size of 5200, the mean NRMSE is 0.148.

Table 1: NARMA10 benchmark results with 400 nodes.  Samples | Mean | Std | Min
---|---|---|---
$2000$ | $0.142$ | $6.47\times 10^{-3}$ | $0.129$
$25000$ | $0.148$ | $5.90\times 10^{-3}$ | $0.139$
Figure 4: NARMA10 benchmark results: (a) shows the amplitude of the target
(red line) and the predicted signal (green line) vs index of 2000 samples and
(b) shows the normalized residue of target and predicted signal.

###  4.2 Santa Fe laser data prediction

The data set A from the 1994 time series prediction competition organized by
the Santa Fe Institute was an time series obtained from measuring a NH3 laser
and it is good example of realistic data [46]. Santa Fe laser data prediction
is one step series prediction and we use 4000 data points in our test case.
The performance is measured using NMSE.

By sweeping the system parameters, the optimum operating point is identified
at gain $G=0.58$, input scaling $\gamma=0.0029$, bias $\phi=0.1\pi$, total of
$N=400$ virtual nodes for a data points of 4000. Results of the experiment is
listed in Table 2. Figure 5 shows plots of the target and prediction, along
with the Residue for the first 1000 samples.

Table 2: Santa Fe laser data single step benchmark results.  Nodes# | Mean | Std | Min
---|---|---|---
$400$ | $2.71\times 10^{-2}$ | $6.20\times 10^{-4}$ | $2.54\times 10^{-2}$
$950$ | $6.73\times 10^{-3}$ | $3.34\times 10^{-4}$ | $5.66\times 10^{-3}$
Figure 5: Santa Fe benchmark results: (a) shows the amplitude of the target
(red line) and the predicted signal (green line) as a function of sample size
and (b) shows the normalized residue of target and prediction.

###  4.3 Isolate spoken digit recognition

Speech recognition is a commonly used benchmark for testing the performance of
neural networks. Acoustic feature are more pronounced in frequency domain
compared to time domain. Hence the input audio data are first pre-processed by
decomposing the time-domain information into frequency-time information. In
our implementation we are using Lyons cochlear ear model to get the
cochleagram which mimics filtering that occurs in nature [47].

Figure 6: Graphical illustration of isolated spoken digit recognition task.
(a) Uniformly distributed input mask with values in range [-1,1] and dimension
$N$x$N_{ch}$ where $N=400$ Nodes and $N_{ch}=77$ (b) Cochleagram generated
from audio file of dimension $N_{ch}$x$N_{T}$ corresponding to digit 9. (c)
Resultant product of input mask and cochleagram is serialized and injected
into reservoir. (d) Output layer weights of dimension $N_{d}$x$N$ where number
of output classes is $N_{d}=10$ (e) The output from reservoir is serially
captured and reshaped to $N$x$N_{T}$, this will give node matrix. (f) Product
of output weight with reservoir state matrix give estimation matrix. Output
class is predicted by getting the maximum argument corresponding to row-wise
mean of estimation matrix.

The input audio data for this experiment is from AudioMNIST database [48].
This is a free and open database that contain 30,000 spoken digits(0-9) audio
samples of 60 different speakers with a sampling rate of 48 kHz. To use these
audio files with our RC, we first down-sample the audio to 12 kHz using the
librosa package [49] and pad audio files with random length zero at beginning
and end to make total sample size of $12,000$. Next we generate the
cochleagram as shown in Fig. 6(b), which is calculated using “lyon 1.0.0"
python package. In this test, we consider first 5 speakers from the AudioMNIST
database to create 20 balanced subset, each containing 5 speaker $\times$ 10
digit $\times$ 2 utterance each, i.e. 100 audio sample. The training and
testing dataset will contains 5 such subset with the first 4 used for training
and the last one for testing. In-order to obtain unbiased result, 50 datasets
are formed using random combinations of 5 subsets from 20 subsets and the word
error rate (WER) is measured for each combinations.

Figure 6 illustrate the flow for spoken digit recognition task. The input mask
is a real valued matrix of dimension $N\times N_{ch}$ with uniform
distribution in [-1,1], where $N$ is number of virtual nodes and $N_{ch}$ is
the number of cochleagram channels. Figure 6(b) shows the cochleagram for
digit 9 with dimension $N_{ch}\times N_{T}$ where $N_{T}$ is the index of new
time representation that depends on decimation factor. As shown in Fig. 6(c),
the resultant product of input mask and cochleagram is serially injected into
reservoir by flattening the matrix. The output of the reservoir is de-
serialized to get the reservoir state matrix of dimensions $N\times N_{T}$, as
shown in Fig. 6(e). The estimation matrix is calculated by multiplying the
output weight matrix with reservoir state matrix. The Output class is
predicted by getting the maximum argument corresponding to row-wise mean of
estimation matrix [50].

After tuning the system parameters and based on their error evaluations, we
found the optimum operating point at gain $G=0.5$, input scaling $\gamma=300$,
bias $\phi=0.35$, total of $N=400$ virtual nodes. Results of experiment are
listed in table 3.

Table 3: WER for Isolated spoken digit recognition task for training and testing phase. | Training WER | Testing WER
---|---|---
Mean | $0.05\%$ | $0.34\%$
Std | $0.13\%$ | $0.51\%$
Table 4: Performance Metric Comparison of Various RC Systems. AWG:Arbitrary
Wavefrom generator, DAQ: Data Acquisition.

Reference | Opto-electronic RC system | Nodes # | Test case | Performance
---|---|---|---|---
[31] |  | Band-pass, optical phase dynamical
---
variable,
AWG @ 24 GS/s.
DAQ @ 80 GS/s
1113 |  | Spoken digit recognition
---
Frequency Channel =86
Training/Testing Size=475/25 digits,
repeated 20 times.
WER=$0.04\%$
[42] |  Low-pass, voltage dynamic variable | 400 |  | Santa Fe laser time-series
---
One step prediction
NMSE=$0.124\pm 4\times 10^{-4}$
400 |  | Spoken digit recognition
---
Frequency Channel=86
Training/Testing Size=475/25 digits,
repeated 20 times
WER=$0.04\%\pm 0.017\%$
[43] |  | Band-pass, voltage dynamic variable,
---
fully analog system,
AWG @ 200 MS/s, 16 bits resolution,
DAQ @ 200 MS/s, 12 bits resolution.
47 |  | NARMA10
---
Training/Testing size=1000/1000.
NMSE=$0.230\pm 0.023$
[51] |  | Band-pass, voltage dynamic variable,
---
Simultaneous task
AWG @ 200 MS/s, 16 bits resolution.
DAQ @ 200 MS/s, 12 bits resolution.
50 |  | NARMA10
---
Training/Testing size=3000/6000.
NMSE=$0.181\pm 0.013$
[52] |  | Band-pass, voltage dynamic variable,
---
optimization via backpropagation
80 |  | NARMA10
---
Back propogation algorithm with
20,000 iterations
NRMSE=$0.185$
[53] | Band-pass, voltage dynamic variable | 400 |  | Santa Fe laser time-series
---
Training/Testing size=3000/1000.
NMSE=$0.02$
[54] |  | Low-pass,
---
wavelength dynamical variable
150 |  | Spoken digit recognition
---
Frequency Channel =86
Training/Testing Size=475/25 digits,
repeated 20 times.
WER=$0.6\%\pm 0.2\%$
[55] |  Band-pass, voltage dynamic variable DAQ @200 MS/s | 50 |  | NARMA10
---
Training/Testing size=1000/1000
NMSE=$0.168\pm 0.015$
200 |  | Spoken digit recognition
---
Frequency Channel=86
Training/Testing Size=475/25 digits,
repeated 5 times.
WER=$0.4\%$
This paper |  Band-pass, Voltage dynamic variable, AWG @ 1 MS/s, 16 bits resolution DAQ @ 1 MS/s, 16 bits resolution. | 400 |  | NARMA10
---
Training/Testing size=1000/1000
NRMSE=$0.142\pm 6.5\times 10^{-3}$
950 |  | Santa Fe laser time-series
---
One step prediction 4000 samples
NMSE=$6.73\times 10^{-3}\pm 3.3\times 10^{-4}$
400 |  | Spoken digit recognition
---
Frequency Channel=77
Training/Testing Size=475/25 digits,
repeated 20 times
WER=$0.34\%\pm 0.51\%$

##  5 Conclusions

We have experimentally demonstrated an opto-electronic RC system using EOM and
FPGA. It takes advantage of electro-optic nonlinear transformation in EOM, and
flexible signal generation and controllable timing, highly programmable signal
filtering, and stable delay logic implementation in FPGA. The resultant system
is a stable and fully functional reservoir computer supporting many nodes,
high connectivity, and easily configurable for multifaceted tasks, while
allowing online training. We have tested it with three benchmark tasks:
NARMA-10, SantaFe laser prediction, and Isolate spoken digit recognition. It
achieves 0.142 NRMSE for NARMA-10, 6.73$\times 10^{-3}$ NMSE for SataFe, and a
WER $\sim$ 0.34% for the speech recognition. These results are compared with
the state of the art in Table 4. As seen, the present RC system tops the
prediction accuracy in the first two test, and is the only reported system to
perform well in all three tasks. Moreover, it is able to precisely predict
25,000 steps in NARMA-10 series with impressive low NRMSE of 0.148. Such high
performance in those various benchmark tests clearly demonstrate the
advantages of our system for complex and versatile tasks.

Our system is currently working at a low speed, which is mostly limited by the
sampling rate of ADC and settling time of DAC. A significant speedup is
expected by choosing a faster ADC and DAC. Also, it is prospective to replace
the existing bulk-optical EOM with photonic integrated chips, where
sophisticated and tailored nonlinear transformations can be realized using
nested optical circuits on a single chip [56, 57]. Finally, an FPGA-based
opto-electronic RC system of this or similar design are easily re-configurable
to accommodate even more neurons, high connectivity, and arbitrary topology.
Also, the performance and robustness of the present RC could be further
improved by online training. Those upgrades will invite important applications
in weather and financial forecasting, real-time information processing, and so
on.

##  6 Disclosures

Disclosures. The authors declare no conflicts of interest.

## References

  * [1] Q. Liu, L. Wang, A. G. Frutos, A. E. Condon, R. M. Corn, and L. M. Smith, “DNA computing on surfaces,” Nature 403, 175–179 (2000).
  * [2] P. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, R. Byer, M. Fejer, H. Mabuchi, and Y. Yamamoto, “A fully programmable 100-spin coherent ising machine with all-to-all connections,” Science 354, 614–617 (2016).
  * [3] R. A. Nawrocki, R. M. Voyles, and S. E. Shaheen, “A mini review of neuromorphic architectures and implementations,” IEEE Transactions on Electron Devices 63, 3819–3829 (2016).
  * [4] F. Amato, A. López, E. M. Peña-Méndez, P. Vaňhara, A. Hampl, and J. Havel, “Artificial neural networks in medical diagnosis,” Journal of Applied Biomedicine 11, 47 – 58 (2013).
  * [5] S. Kumar, H. Zhang, and Y.-P. Huang, “Large-scale ising emulation with four body interaction and all-to-all connections,” Communications Physics 3, 108 (2020).
  * [6] A. Krizhevsky, I. Sutskever, and G. Hinton, “Imagenet classification with deep convolutional neural networks,” Neural Information Processing Systems 25 (2012).
  * [7] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” Nature 521, 436–444 (2015).
  * [8] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, S. Petersen, C. Beattie, A. Sadik, I. Antonoglou, H. King, D. Kumaran, D. Wierstra, S. Legg, and D. Hassabis, “Human-level control through deep reinforcement learning,” Nature 518, 529–533 (2015).
  * [9] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. van den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S. Dieleman, D. Grewe, J. Nham, N. Kalchbrenner, I. Sutskever, T. Lillicrap, M. Leach, K. Kavukcuoglu, T. Graepel, and D. Hassabis, “Mastering the game of go with deep neural networks and tree search,” Nature 529, 484–489 (2016).
  * [10] Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, and M. Soljačić, “Deep learning with coherent nanophotonic circuits,” Nature Photonics 11, 441 (2017).
  * [11] J. Feldmann, N. Youngblood, C. D. Wright, H. Bhaskaran, and P. W. H. P., “All-optical spiking neurosynaptic networks with self-learning capabilities,” Nature 569, 208 (2019).
  * [12] X. Sui, Q. Wu, J. Liu, Q. Chen, and G. Gu, “A review of optical neural networks,” IEEE Access 8, 70773–70783 (2020).
  * [13] T. Bu, S. Kumar, H. Zhang, I. Huang, and Y.-P. Huang, “Single-pixel pattern recognition with coherent nonlinear optics,” Optics Letters 45, 6771–6774 (2020).
  * [14] Y. Zuo, B. Li, Y. Zhao, Y. Jiang, Y.-C. Chen, P. Chen, G.-B. Jo, J. Liu, and S. Du, “All-optical neural network with nonlinear activation functions,” Optica 6, 1132–1137 (2019).
  * [15] T. Zhou, L. Fang, T. Yan, J. Wu, Y. Li, J. Fan, H. Wu, X. Lin, and Q. Dai, “In situ optical backpropagation training of diffractive optical neural networks,” Photon. Res. 8, 940–953 (2020).
  * [16] D. Verstraeten, B. Schrauwen, M. D’Haene, and D. Stroobandt, “An experimental unification of reservoir computing methods,” Neural networks : the official journal of the International Neural Network Society 20, 391—403 (2007).
  * [17] L. Appeltant, M. C. Soriano, G. Van der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. R. Mirasso, and I. Fischer, “Information processing using a single dynamical node as complex system,” Nature Communications 2, 468 (2011).
  * [18] P. Antonik, M. Haelterman, and S. Massar, “Brain-Inspired Photonic Signal Processor for Generating Periodic Patterns and Emulating Chaotic Systems,” Physical Review Applied 7, 054014 (2017).
  * [19] D. Brunner, B. Penkovsky, B. A. Marquez, M. Jacquot, I. Fischer, and L. Larger, “Tutorial: Photonic neural networks in delay systems,” Journal of Applied Physics 124, 152004 (2018).
  * [20] Y. K. Chembo, “Machine learning based on reservoir computing with time-delayed optoelectronic and photonic systems,” Chaos 30, 013111 (2020).
  * [21] W. Maass, T. Natschläger, and H. Markram, “Real-time computing without stable states: A new framework for neural computation based on perturbations,” Neural Computation 14, 2531–2560 (2002).
  * [22] H. Jaeger and H. Haas, “Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication,” Science 304, 78–80 (2004).
  * [23] K. Vandoorne, W. Dierckx, B. Schrauwen, D. Verstraeten, R. Baets, P. Bienstman, and J. Van Campenhout, “Toward optical signal processing using Photonic Reservoir Computing,” Optics Express 16, 11182 (2008).
  * [24] K. Vandoorne, J. Dambre, D. Verstraeten, B. Schrauwen, and P. Bienstman, “Parallel Reservoir Computing Using Optical Amplifiers,” IEEE Transactions on Neural Networks 22, 1469–1481 (2011).
  * [25] F. Duport, B. Schneider, A. Smerieri, M. Haelterman, and S. Massar, “All-optical reservoir computing,” Optics Express 20, 22783 (2012).
  * [26] D. Brunner, M. C. Soriano, C. R. Mirasso, and I. Fischer, “Parallel photonic information processing at gigabyte per second data rates using transient states,” Nature Communications 4, 1364 (2013).
  * [27] C. Mesaritakis, V. Papataxiarhis, and D. Syvridis, “Micro ring resonators as building blocks for an all-optical high-speed reservoir-computing bit-pattern-recognition system,” Journal of the Optical Society of America B 30, 3048 (2013).
  * [28] A. Dejonckheere, F. Duport, A. Smerieri, L. Fang, J.-L. Oudar, M. Haelterman, and S. Massar, “All-optical reservoir computer based on saturation of absorption,” Optics Express 22, 10868 (2014).
  * [29] T. Yamane, Y. Katayama, R. Nakane, G. Tanaka, and D. Nakano, “Wave-Based Reservoir Computing by Synchronization of Coupled Oscillators,” in _Neural Information Processing,_ vol. 9491 S. Arik, T. Huang, W. K. Lai, and Q. Liu, eds. (Springer International Publishing, Cham, 2015), pp. 198–205. Series Title: Lecture Notes in Computer Science.
  * [30] C. Mesaritakis, A. Bogris, A. Kapsalis, and D. Syvridis, “High-speed all-optical pattern recognition of dispersive Fourier images through a photonic reservoir computing subsystem,” Optics Letters 40, 3416 (2015).
  * [31] L. Larger, A. Baylón-Fuentes, R. Martinenghi, V. S. Udaltsov, Y. K. Chembo, and M. Jacquot, “High-Speed Photonic Reservoir Computing Using a Time-Delay-Based Architecture: Million Words per Second Classification,” Physical Review X 7, 011015 (2017).
  * [32] G. Tanaka, R. Nakane, T. Yamane, S. Takeda, D. Nakano, S. Nakagawa, and A. Hirose, “Waveform classification by memristive reservoir computing,” in _Neural Information Processing,_ vol. 10637 D. Liu, S. Xie, Y. Li, D. Zhao, and E.-S. M. El-Alfy, eds. (Springer International Publishing, Cham, 2017), pp. 457–465. Series Title: Lecture Notes in Computer Science.
  * [33] G. Van der Sande, D. Brunner, and M. C. Soriano, “Advances in photonic reservoir computing,” Nanophotonics 6, 561–576 (2017).
  * [34] Y. Hou, G. Xia, W. Yang, D. Wang, E. Jayaprasath, Z. Jiang, C. Hu, and Z. Wu, “Prediction performance of reservoir computing system based on a semiconductor laser subject to double optical feedback and optical injection,” Optics Express 26, 10211 (2018).
  * [35] A. Lugnan, A. Katumba, F. Laporte, M. Freiberger, S. Sackesyn, C. Ma, E. Gooskens, J. Dambre, and P. Bienstman, “Photonic neuromorphic information processing and reservoir computing,” APL Photonics 5, 020901 (2020).
  * [36] D. Canaday, A. Griffith, and D. J. Gauthier, “Rapid time series prediction with a hardware-based reservoir computer,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 123119 (2018).
  * [37] B. Penkovsky, L. Larger, and D. Brunner, “Efficient design of hardware-enabled reservoir computing in fpgas,” Journal of Applied Physics 124, 162101 (2018).
  * [38] I. Estébanez, I. Fischer, and M. C. Soriano, “Constructive role of noise for high-quality replication of chaotic attractor dynamics using a hardware-based reservoir computer,” Phys. Rev. Applied 12, 034058 (2019).
  * [39] F. Böhm, G. Verschaffelt, and G. Van der Sande, “A poor man’s coherent ising machine based on opto-electronic feedback systems for solving optimization problems,” Nature Communications 10, 3538 (2019).
  * [40] J. D. Hart, D. C. Schmadel, T. E. Murphy, and R. Roy, “Experiments with arbitrary networks in time-multiplexed delay systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 121103 (2017).
  * [41] J. D. Hart, R. Roy, D. Müller-Bender, A. Otto, and G. Radons, “Laminar chaos in experiments: Nonlinear systems with time-varying delays and noise,” Phys. Rev. Lett. 123, 154101 (2019).
  * [42] L. Larger, M. C. Soriano, D. Brunner, L. Appeltant, J. M. Gutierrez, L. Pesquera, C. R. Mirasso, and I. Fischer, “Photonic information processing beyond Turing: an optoelectronic implementation of reservoir computing,” Optics Express 20, 3241 (2012).
  * [43] F. Duport, A. Smerieri, A. Akrout, M. Haelterman, and S. Massar, “Fully analogue photonic reservoir computer,” Scientific Reports 6, 22381 (2016).
  * [44] X. Bao, Q. Zhao, and H. Yin, “Efficient optoelectronic reservoir computing with three-route input based on optical delay lines,” Applied Optics 58, 4111 (2019).
  * [45] H. Jaeger, “Adaptive nonlinear system identification with echo state networks,” in _NIPS,_ (2002).
  * [46] C. Chatfield and A. S. Weigend, “Time series prediction: Forecasting the future and understanding the past: Neil a. gershenfeld and andreas s. weigend, 1994, ’the future of time series’, in: A.s. weigend and n.a. gershenfeld, eds., (addison-wesley, reading, ma), 1-70,” International Journal of Forecasting 10, 161–163 (1994).
  * [47] R. Lyon, “A computational model of filtering, detection, and compression in the cochlea,” in _ICASSP ’82. IEEE International Conference on Acoustics, Speech, and Signal Processing,_ vol. 7 (1982), pp. 1282–1285.
  * [48] S. Becker, M. Ackermann, S. Lapuschkin, K.-R. Müller, and W. Samek, “Interpreting and explaining deep neural networks for classification of audio signals,” CoRR abs/1807.03418 (2018).
  * [49] B. McFee, V. Lostanlen, A. Metsai, M. McVicar, S. Balke, C. Thomé, C. Raffel, F. Zalkow, A. Malek, Dana, K. Lee, O. Nieto, J. Mason, D. Ellis, E. Battenberg, S. Seyfarth, R. Yamamoto, K. Choi, viktorandreevichmorozov, J. Moore, R. Bittner, S. Hidaka, Z. Wei, nullmightybofo, D. Hereñú, F.-R. Stöter, P. Friesch, A. Weiss, M. Vollrath, and T. Kim, “librosa/librosa: 0.8.0,” (2020).
  * [50] F. Abreu Araujo, M. Riou, J. Torrejon, S. Tsunegi, D. Querlioz, K. Yakushiji, A. Fukushima, H. Kubota, S. Yuasa, M. D. Stiles, and J. Grollier, “Role of non-linear data processing on speech recognition task in the framework of reservoir computing,” Scientific Reports 10, 328 (2020).
  * [51] F. Duport, A. Smerieri, A. Akrout, M. Haelterman, and S. Massar, “Virtualization of a photonic reservoir computer,” Journal of Lightwave Technology 34, 2085–2091 (2016).
  * [52] M. Hermans, P. Antonik, M. Haelterman, and S. Massar, “Embodiment of learning in electro-optical signal processors,” Phys. Rev. Lett. 117, 128301 (2016).
  * [53] M. C. Soriano, S. Ortín, D. Brunner, L. Larger, C. R. Mirasso, I. Fischer, and L. Pesquera, “Optoelectronic reservoir computing: tackling noise-induced performance degradation,” Optics Express 21, 12 (2013).
  * [54] R. Martinenghi, S. Rybalko, M. Jacquot, Y. K. Chembo, and L. Larger, “Photonic Nonlinear Transient Computing with Multiple-Delay Wavelength Dynamics,” Physical Review Letters 108, 244101 (2012).
  * [55] Y. Paquot, F. Duport, A. Smerieri, J. Dambre, B. Schrauwen, M. Haelterman, and S. Massar, “Optoelectronic Reservoir Computing,” Scientific Reports 2, 287 (2012).
  * [56] J.-Y. Chen, Z.-H. Ma, Y. M. Sua, Z. Li, C. Tang, and Y.-P. Huang, “Ultra-efficient frequency conversion in quasi-phase-matched lithium niobate microrings,” Optica 6, 1244 (2019).
  * [57] M. Jin, J.-Y. Chen, Y. M. Sua, and Y.-P. Huang, “High-extinction electro-optic modulation on lithium niobate thin film,” Opt. Lett. 44, 1265–1268 (2019).