# Online Training of Spiking Recurrent Neural Networks with Phase-Change
Memory Synapses

Yiğit Demirağ
The Institute of Neuroinformatics
University of Zurich, ETH Zurich
Zürich, Switzerland
<EMAIL_ADDRESS>
&Charlotte Frenkel
The Institute of Neuroinformatics
University of Zurich, ETH Zurich
Zürich, Switzerland
&Melika Payvand
The Institute of Neuroinformatics
University of Zurich, ETH Zurich
Zürich, Switzerland
&Giacomo Indiveri
The Institute of Neuroinformatics
University of Zurich, ETH Zurich
Zürich, Switzerland

###### Abstract

Spiking Recurrent Neural Networks are a promising tool for solving a wide
variety of complex cognitive and motor tasks, due to their rich temporal
dynamics and sparse processing. However training spiking RNNs on dedicated
neuromorphic hardware is still an open challenge. This is due mainly to the
lack of local, hardware-friendly learning mechanisms that can solve the
temporal credit assignment problem and ensure stable network dynamics, even
when the weight resolution is limited. These challenges are further
accentuated, if one resorts to using memristive devices for in-memory
computing to resolve the von-Neumann bottleneck problem, at the expense of a
substantial increase in variability in both the computation and the working
memory of the spiking RNNs. To address these challenges and enable online
learning in memristive neuromorphic RNNs, we present a simulation framework of
differential-architecture crossbar arrays based on an accurate and
comprehensive Phase Change Memory (PCM) device model. We train a spiking RNN
whose weights are emulated in the presented simulation framework, using a
recently proposed e-prop learning rule. Although e-prop locally approximates
the ideal synaptic updates, it is difficult to implement the updates on the
memristive substrate due to substantial PCM non-idealities. We compare several
widely adapted weight update schemes that primarily aim to cope with these
device non-idealities and demonstrate that accumulating gradients can enable
online and efficient training of spiking RNNs on memristive substrates.

##  1 Introduction

RNNs are exceptionally expressive Khrulkov_etal17  class of neural networks
that have been successfully adapted in many domains such as audio processing,
optical flow, language modeling and Reinforcement Learning (RL) Oord_etal16 ;
Teed_Deng20 ; Brown_etal20 ; Berner_etal19 ; Ha_Schmidhuber18 ; Vinyals_etal19
. The power of RNNs lie in their architecture that allows processing of long
and complex sequential data. Due to its recurrent architecture, each neuron
contribute to the network processing at various points of the computation. The
resulting efficiency of RNNs is also evident in the mammalian brain with
massive lateral recurrent connections in neocortex Douglas_Martin98 ;
Douglas_etal95a . However, the training of RNN topologies is notably difficult
under constrained memory and computational resources Zenke_Neftci21 .

Current hardware implementations of neural networks still cannot compete with
the energy efficiency of biological systems. One of the main reasons for this
difference is due to the data movement between the separated processing and
memory units of the von-Neumann architectures used to implement the artificial
neural networks. Recently, new types of compact nanoscale devices have
garnered significant attention for implementing artificial synapses
Ambrogio_etal18 ; Li_etal18a ; Prezioso_etal15 ; dalgaty_etal21 ; cai_etal20 ;
Sebastian_etal19  that can implement in-situ learning and break the von
Neumann bottleneck Backus78 ; Indiveri_Liu15 . These memristive devices are
particularly promising for use in Spiking Neural Network (SNN) architectures,
especially for low-power, sparse and massively parallel spike-based
neuromorphic systems, that emulate the principles of computation observed in
biological brains Payvand_etal19 ; Chicca_Indiveri20 . In these systems, the
synapses (memory) and neurons (processing units) are arranged in a crossbar
architecture as shown in Fig. 1a, where memristive devices are placed at the
junctions, to store the synaptic efficacy in their programmable conductance
values. The crossbar architecture implements in-memory computation of synaptic
propagation between neurons by a single physics-based operation following
Ohm’s Law and Kirchhoff’s Law, intrinsically supporting the sparse and event-
driven nature of SNNs. As recently shown in Peng_etal19  and Peng_etal20  for
32 nm technology, compared to a digital Complementary Metal-Oxide-
Semiconductor (CMOS) implementation (i.e. Static Random Access Memory (SRAM))
such memristive crossbar arrays enable a denser solution, and have
respectively lower and similar dynamic energy consumption during inference and
training. Moreover, the non-volatile nature of memristive devices reduces the
static power consumption usually linked with volatile CMOS memory storage.
Therefore, in-memory acceleration of spiking RNNs with non-volatile multi-bit-
resolution memristive devices is a promising direction for scalable
neuromorphic hardware for temporal signal processing.

PCM devices are among the most mature emerging resistive memory technologies.
Their tiny footprint, fast read/write operation and multi-bit storage capacity
make PCMs an ideal candidate for implementing in-memory computation of
synaptic propagation Burr_etal16a ; Burr_etal17 . Consequently, there has been
an increased interest for the employment of PCM technology in neuromorphic
computing applications Ambrogio_etal18 ; Tuma_etal16 ; Karunaratne_etal20a ;
Demirag_etal21a .

Single PCM device can have 3-4 bits of resolution Gallo_etal18 , however,
unlike conventional digital CMOS memories, they are affected by severe non-
idealities as their switching physics is governed by inherently stochastic
Joule heating. Moreover, molecular dynamics of PCMs give rise to the $1/f$
noise behavior and the structural relaxation, resulting in cycle-to-cycle
variation, in addition to their device-to-device variability present due to
fabrication effects.

To account for these device non-idealities, one could follow a hardware-
algorithm co-design approach with chip-in-the-loop setups that leverage PCM
crossbar hardware Ambrogio_etal18 . However, training a neural network
requires the iterative evaluation of multiple network architectures,
modifications to the learning rule and tuning hyperparameters on large
datasets, which are extremely time/resource consuming with chip-in-the-loop
setups. On the other hand, a software network simulation framework with a
statistical model of memristive devices offers much faster iteration times and
a better understanding of device effects on the training, due to increased
observability of internal state variables. However, it is extremely important
to have a very accurate statistical model of the devices simulated, to
optimize the training procedure of the network before moving to chip-in-the-
loop training.

In this paper, we investigate whether a spiking RNN can be trained with a
local learning rule despite the adverse impacts of in-memory computing with
memristive devices such as write noise, read noise, temporal conductance
evolution (i.e., drift) and the limited bit precision. To do so, we first
build on the statistical PCM model from Nandakumar et al. Nandakumar_etal18
to faithfully model a crossbar array comprising a differential memristor
configuration (Section 2.1). Then, we define a target spiking RNN architecture
and describe the properties of an ideal learning rule and select e-prop
algorithm Bellec_etal20  to train the network (Section 2.2). In order to map
ideal weight updates calculated by e-prop to memristor conductances on the
crossbar array, we implement multiple memristor-aware weight update methods
that are optimized to cope with device non-idealities (Section 2.3). Finally,
we report a training scheme for spiking RNNs which exploits in-memory
computing with extremely sparse activity and a reduced number of conductance
updates for energy efficient training (Section 4).

##  2 Building Blocks for In-Memory Acceleration of RNN Training for
Neuromorphic Processors

In the following, we describe the main components of our simulation framework
for training spiking RNNs with PCM synapses 111The code is available at
https://github.com/YigitDemirag/srnn-pcm.

###  2.1 PCM Synapses

The material design of a nanoscale PCM device typically includes the switching
material $\text{Ge}_{2}\text{Sb}_{2}\text{Te}_{5}$ (GST) placed between two
metal electrode layers forming a mushroom-type structure (Fig 1b). By applying
short electrical pulses, the temperature distribution inside the PCM device
can be momentarily modified via Joule heating. The controlled temperature
levels can switch the molecular configuration of GST between the amorphous
(high-resistance) and crystalline (low-resistance) states Demirag18 .

To increase the volume of the amorphous state in the device, a short and high-
amplitude electrical pulse (RESET pulse) is applied to the device terminals.
Increasing the temperature to around 900 K melts a significant portion of the
GST, and if the GST quenches rapidly, melted region forms an amorphous
configuration. On the contrary, to increase the crystalline volume, typically
a longer and smaller-amplitude electrical pulse (SET pulse) is applied to the
device terminals. In this case, temperature rises to around 400-600 K, which
initiates the growth of available crystal nuclei inside the GST and increases
the crystalline volume. To read the device conductance, an electrical pulse
with an even smaller amplitude (READ pulse) is applied to the device, in order
to prevent any phase transition. The amplitude and duration of SET, RESET and
READ pulses depend on the GST composition and device volume. Nevertheless for
a typical mushroom type geometry with 100 nm GST and 20 nm heater radius (see
Fig. 1b), a SET pulse of 2.5V amplitude and 100 ns duration, a RESET pulse of
3.5 V amplitude and 20 ns duration and READ pulse of 0.2 V 50 ns duration can
be used Demirag18 .

Figure 1: a. PCM devices can be arranged in a crossbar architecture to emulate
both a non-volatile synaptic memory and a parallel and asynchronous synaptic
propagation using in-memory computation. b. Mushroom-type geometry of a single
PCM device. The conductance of the device can be reconfigured by changing the
volume ratio of amorphous and crystalline regions.

In practice however, the programming operations of such devices suffer from
write noise, read noise, and the electrical conductance drift
Gallo_Sebastian20 . The asymmetry of the SET and RESET operations, and the
nonlinear conductance response with respect to the number and frequency of the
pulses applied further complicates the precise programming of the device
conductance to target values. It is therefore critical to fully capture these
non-ideal phenomena and device dynamics in the network models, because only
then important metrics such as the robustness of the weight update rule,
hyperparameter choices or the training duration can be realistically
evaluated, reflecting the real-world deployment of the neuromorphic system.
Many comprehensive models have been proposed to describe electrical
Gallo_etal15 , thermal Gallo_etal16 , structural Ielmini_etal07 ;
Karpov_etal07  and phase-change Redaelli_etal08 ; Salinga_etal13  properties
of PCMs. These device models either require solving on-the-fly differential
equations whose numerical convergence is not guaranteed, or do not incorporate
inter- and intra- device stochasticity, or are designed for pulse shapes and
current-voltage sweeps that do not reflect the operational conditions on the
circuit Demirag18 .

In this study, we selected the statistical device model from Nandakumar et al.
Nandakumar_etal18 , which comprises all major forms of PCM-specific non-
idealities based on measurements from 10,000 devices. The model includes the
nonlinear conductance change with respect to applied pulses, conductance-
dependent write and read operation stochasticity, and the temporal drift
effect (Fig. 2). The model keeps a programming history variable, which
represents the nonlinear device response to consecutive SET pulses and is
updated after the application of each pulse. Following the application of a
new SET pulse, the model samples the conductance change $\Delta G$ from a
Gaussian distribution whose mean and standard deviation is based on the
programming history and the previous conductance. The drift behavior is then
included following the empirical exponential drift model Karpov_etal07
$G(t)=G\left(T_{0}\right)\left(t/T_{0}\right)^{-v}$, where $G(T_{0})$ is the
estimated conductance after a WRITE pulse at time $T_{0}$ and $G(t)$ is the
final device conductance considering the effect of the drift. Additionally,
the model takes into account the $1/f$ READ noise Nardone_etal09 , which
increases monotonically with the device conductance. Overall, this statistical
PCM model captures the stochastic conductance changes due to the application
of SET and READ pulses and estimates the temporal conductance drift arising
from the structural relaxation.

In order to integrate this model into neural network simulations, we developed
a comprehensive PCM crossbar array simulation framework in PyTorch
Paszke_etal19a . Our crossbar array simulation framework can keep track of all
simulated PCM devices simultaneously, enabling realistic SET, RESET and READ
operations (for implementation details, see Supplementary Note 1). Section 3
will describe how this framework is used to represent synaptic weights of a
RNN.

Figure 2: The chosen PCM model from Nandakumar_etal18  captures the major
device non-idealities. a. The WRITE model enables calculation of the
conductance increase with each consecutive SET pulse applied to the device.
The band illustrates one standard deviation. b. The READ model enables
calculation of $1/f$ noise, increasing as a function of conductance. c. The
DRIFT model calculates the temporal conductance evolution as a function time.
$T_{0}$ indicates the time of measurement after the initial programming of the
device.

###  2.2 Credit Assignment Solutions for Recurrent Network Architectures

The credit assignment problem refers to the problem of determining the amount
of change required in each synaptic weight to achieve the network’s desired
behavior Richards_etal19 . The nature of this problem is intertwined with
network architecture which describes the arrangement of synapses, neurons and
their operational principles. Hence many proposed solutions to the credit
assignment problem in SNNs landscape are specific to a network architecture
components e.g., eligibility traces Fremaux_Gerstner16 , dendritic
Sacramento_etal18  or neuromodulatory Pozzi_etal18  signals.

In this work, we simulated the neuron dynamics using the Leaky Integrate-and-
Fire (LIF) neuron model Gerstner_Kistler02 ; Dayan_Abbott01 . LIF neurons are
stateful through their membrane potential, and their temporal dynamics in a
spiking RNN can be simulated with the following discrete-time equations
Bellec_etal20 :

$\displaystyle v_{j}^{t+1}$ $\displaystyle=\alpha v_{j}^{t}+\sum_{i\neq
j}W_{ji}^{\mathrm{rec}}z_{i}^{t}+\sum_{i}W_{ji}^{\mathrm{in}}x_{i}^{t}-z_{j}^{t}v_{\mathrm{th}}$
(1) $\displaystyle z_{j}^{t}$
$\displaystyle=H\left(\frac{v_{j}^{t}-v_{\mathrm{th}}}{v_{\mathrm{th}}}\right)$
$\displaystyle y_{k}^{t+1}$ $\displaystyle=\kappa
y_{k}^{t}+\sum_{j}W_{kj}^{\mathrm{out}}z_{j}^{t}$

where $v_{j}^{t}$ is the membrane voltage of neuron $j$ at time $t$. The
output state of a neuron is a binary variable, $z_{j}^{t}$, that can either
indicate a spike, $1$, or no spike, $0$. The neuron spikes when the membrane
voltage exceeds the threshold voltage $v_{th}$, a condition that is
implemented based on the Heaviside function $H$. The parameter
$\alpha\in[0,1]$ is the membrane decay factor calculated as $\alpha=e^{-\delta
t/\tau_{m}}$, where $\delta t$ is the discrete time step resolution of the
simulation and $\tau_{m}$ is the neuronal membrane decay time constant,
typically tens of milliseconds. The network activity is driven by input spikes
$x_{i}^{t}$. Input and recurrent weights are represented as
$W_{ji}^{\mathrm{in}}$ and $W_{ji}^{\mathrm{rec}}$ respectively. At the output
layer, the recurrent spikes are fed through readout weights
$W_{kj}^{\mathrm{out}}$ into a single layer of leaky integrator units $y_{k}$
with the decay factor $\kappa\in[0,1]$. This continuous valued output unit is
analogous to a motor function which generates coherent motor output patterns
Sussillo_Abbott09  of the type shown in Fig. 3.

The aim of the training process is to find optimal network weights
$\\{W_{ji}^{\mathrm{inp}},W_{ji}^{\mathrm{rec}}$ and
$W_{kj}^{\mathrm{out}}\\}$, that maximize the performance of the network
Richards_etal19 . One approach is to train the spiking RNNs offline, using any
optimization method, and then transfer the weights to the crossbar via an
iterative procedure for inference applications Diehl_etal15 . However this
would not take into account the non-idealities of PCM synapses. Our aim is to
develop an online training framework that can compensate for these non-
idealities. To be suitable for neuromorphic hardware, this framework must use
a learning algorithm that is (i) local, (ii) online and (iii) well-tested
beyond toy problems. For example, the FORCE algorithm performs well on motor
tasks Nicola_Clopath17a ; Sussillo_Abbott09 , however the weight updates
require the knowledge of all synaptic weights, which violates the first
requirement. Similarly, training the network with Back-propagation Through
Time (BPTT) using surrogate gradients Neftci_etal19 ; Lee_etal16  is not an
option as it requires buffering all intermediate neuron states to calculate
updates backward in time, which violates the second requirement.

One promising solution to the credit assignment problem in the spiking RNN
landscape is provided by the e-prop algorithm Bellec_etal20 . It has been
shown that e-prop can offer an accuracy similar to that of Long Short Term
Memory (LSTM) networks Hochreiter_Schmidhuber97  trained with BPTT on complex
temporal tasks. As this algorithm works by factorizing the gradients of BPTT
as a sum of products between instantaneous learning signals and local
eligibility traces, it is both online and local. Similar to e-prop, the
recently published Online Spatio-Temporal Learning (OSTL) algorithm
Bohnstingl_etal20  also supports complex recurrent architectures, and in
addition it is able to generate the exact gradients calculated with BPTT.

###  2.3 Memristor-Aware Weight Update Optimization

Many neuromorphic processors with on-chip learning include a Learning Block
(LB) connected to the neuron circuits Qiao_etal15 ; Payvand_etal18 ;
Payvand_etal19 . The function of LB is to continuously compare the activity of
the neuron with a desired target activity, and (ii) estimate the amount of
required change in the neuron’s synaptic weights, based on the implemented
learning rule and error function. When the weight update condition is met,
which can be on the arrival of the pre-synaptic spike, on the post-synaptic
spike, or on a error-driven signal Payvand_etal20 , the LB instructs the
corresponding synapses to update the weight with the calculated weight change.
However, memristive synapses cannot be easily programmed to target specific
conductances, due to device non-idealities such as limited bit precision,
asymmetric SET/RESET updates and programming noise.

Weight update mechanisms in memristive architectures are primarily designed to
cope with memristive non-idealities to execute the desired weight changes on
synapses. These mechanisms are typically single-shot, i.e. one or multiple
gradual SET pulses are applied without reading the device state. This avoids
iterative write-read-verify schemes or modifications of the pulse shape during
the training, and hence enables simpler circuits with a better energy
efficiency. In our framework we use a differential synaptic configuration
Nandakumar_etal18 ; Boybat_etal18 ; Nandakumar_etal20a , where every synapse
keeps two sets of memristors ($G^{+}$ and $G^{-}$) whose difference gives the
effective synaptic conductance (Fig. 1a), so that both positive and negative
synaptic weights can be represented using memristor conductances.

With this differential configurations however the unidirectional updates used
to update the synaptic weight may result in the saturation of either or both
of the devices in the synapse Nandakumar_etal18 ; Boybat_etal18 ;
Nandakumar_etal20a . One solution to the problem would be to employ a push-
pull mechanism: to increase the synaptic efficacy, the positive memristor
conductance is increased while the negative memristor conductance is decreased
Nair_etal17 . Unfortunately this mechanism is not compatible with PCM devices
as the melt-quenching-based RESET is an abrupt process Burr_etal14 , thus
requiring a refresh mechanism by resetting both positive and negative
memristors and reprogramming them to the their effective conductances when
specific criteria is met. In the following, we will describe four weight
update methods that the neuromorphic processors have widely adopted and are
implemented in the PCM crossbar array simulation framework.

#### The Sign Gradient Descent.

In this method, synaptic weights are updated according to the sign of the
gradient of the loss function, estimated by the learning algorithm. The idea
of the sign gradient descent is to take a fixed step ($\delta$) in the
direction of descending gradient, but neglecting its magnitude. Nevertheless,
under some assumptions, the convergence is guaranteed Balles_etal20a . The
synaptic weight updates $\Delta W$ can be computed as

$\Delta W=-\delta\operatorname{sign}(\nabla\mathcal{L})SL,$ (2)

where $\delta$ is the amount of change, $\nabla\mathcal{L}$ is the gradient
calculated by the chosen learning rule, and $SL=|\nabla{\mathcal{L}}|>\theta$
is a binary variable indicating a stop-learning regime. The latter increases
the stability during learning, which enables updates only when the magnitude
of the gradient is higher than a fixed threshold $\theta$.

Due to its simplicity, the sign gradient descent is popular among memristive
neuromorphic systems Nair_Dudek15 ; Muller_etal17 ; Payvand_etal20a . When
implementing it with memristor synapses, the LB sends a single UP (or DOWN)
pulse to instruct an increase (or decrease) of the synaptic weights. Hence, on
the onset of the weight update, a single SET pulse is applied either to the
$G^{+}$ or the $G^{-}$ PCM device, determined by the sign of the gradient.
However, the effective value of $\delta$ is not a constant due to the WRITE
noise, and not symmetric because SET operation in PCM is gradual whereas RESET
is abrupt.

#### Stochastic Update.

Conventional optimization methods demand the amount of weight change at every
weight update to be 3-4 orders of magnitude smaller than the original weight
Kingma_Ba14 , usually requiring at least an 8-bit weight representation,
whereas a typical PCM device can only represent 3-4 bits of information
Athmanathan_etal16a .

Therefore, even a single WRITE pulse applied to a PCM device leads to an
order-of-magnitude overshoot compared to the desired amount of weight change
typically calculated by the LB. One method used by Nandakumar et al.
Nandakumar_etal18  aims at overcoming this scaling disparity by stochastically
executing (or skipping) the weight updates depending on the magnitude of the
gradient.

We implemented the stochastic weight update on our PCM crossbar simulation by
scaling the gradient magnitude with a scaling factor, $p$, to represent the
update probability such that

$P(\text{update})=\frac{|\nabla\mathcal{L}|}{p}$ (3)

By tuning the scaling factor, the number of devices that get programmed can be
controlled. Single pulses are applied accordingly to the synapses. The refresh
criteria, unlike the original work Nandakumar_etal18 , is checked and applied
before the weight update to prevent the execution of the update on saturated
devices.

#### Multi-memristor Update.

To increase the dynamic range and the resolution of the synaptic weights, the
concept of multiple PCM devices per synapses has been proposed in
Boybat_etal18 , where each synapse consists of $2N$ memristors ($N$ positive
and negative pairs in a differential configuration), the total conductance
determining the synaptic efficacy. The updates are applied to the devices
sequentially, which in turn decreases the smallest mean weight change by a
factor of $2N$ and reduces the variance due to WRITE noise by a factor of
$\sqrt{2N}$ Boybat_etal18  (see Supplementary Note 3).

Our multi-memristor implementation first estimates the the number of pulses
required to match the desired conductance change, by assuming each SET pulse
increases the conductance linearly by $0.75$ $\mu$S as the model outlined in
Section 2.1 suggests (a more realistic approximation of the number of pulses
is available in Nandakumar_etal18 ). Then these SET pulses are applied
sequentially among $N$ PCM devices in a circular queue. The refresh is applied
if any of the memristor pair conductances is higher than 9 $\mu$S and their
difference is less than $4.5\mu$S.

#### Mixed-precision Update.

One solution to close the gap between the update resolution requested by the
LB and the minimum programmable conductance change of PCMs is to accumulate
the updates on a high-precision co-processor until they become reliably
transferable to PCMs Nandakumar_etal20a . Notably, this scheme also
corresponds to the quantization-aware training techniques conventionally used
for the training of quantized artificial neural networks Hubara_etal16 .

In our simulations, the LB accumulates the gradients in a double-precision
floating memory until they are an integer multiple of PCM update granularity
(corresponding to $0.75$ $\mu$S). Then they are converted to a number of
pulses and applied to the memristors. The refresh is applied if either one of
the memristor pair conductance is higher than 9 $\mu$S and their difference is
less than $4.5\mu$S.

##  3 Training a Spiking RNN on a PCM Crossbar Simulation Framework

Figure 3: Overview of the spiking RNN training framework of the proposed PCM
crossbar array simulation framework, illustrated here for a pattern generation
task. The network weights are allocated from three crossbar array models,
$G_{inp}$, $G_{rec}$, $G_{out}$. The network-generated pattern and the target
pattern are compared to produce the learning signal which is fed back to every
neuron. The LB calculates instantaneous weight changes $\Delta W$ using the
e-prop learning rule and has to efficiently transfer the desired weight change
to a conductance change, i.e. $\Delta W\rightarrow\Delta G$, considering PCM
non-idealities.

We used the model of the PCM crossbar array to determine realistic values of
the network parameters $\\{W_{ji}^{\mathrm{in}}$, $W_{ji}^{\mathrm{rec}}$ and
$W_{kj}^{\mathrm{out}}\\}$. In order to represent synaptic weights,
$W\in[-1,1]$ , with the conductance values of PCM devices, $G\in[0.1,12]$
$\mu$S Nandakumar_etal18 , we used the linear relationship
$W=\beta[\sum_{N}G^{+}-\sum_{N}G^{-}]$, where $\sum_{N}G^{+}$ and
$\sum_{N}G^{-}$ are the total conductance of $N$ memristors222$N=1$ for all
weight update methods, except when the multi-memristor update is being used.
representing the potentiation and the depression of the synapse respectively
Nandakumar_etal18 . At this stage, the forward computation (inference) of Eq.
1 is simulated using the PCM crossbar simulation framework that includes the
effects of READ noise and temporal conductance drift.

The weight updates calculated by the e-prop algorithm are applied to the PCM-
based crossbar arrays using each of the methods described in Section 2.3.

##  4 Results

We validated the network using a classical one-dimensional continuous pattern
generation task relevant for many domains including motor control or value
function estimation in RLSussillo_Abbott09  The training dataset is the
slightly modified version of the one-second-long patterns described in
Bellec_etal20  (see Methods for task and training details).

Table 1 summarizes the training performances of spiking RNNs utilizing
different weight update methods on PCM crossbar arrays to realize the target
weight change calculated by the e-prop algorithm. Out of five configurations,
only the mixed-precision approach resulted in an acceptable performance on the
pattern generation task (MSE loss < 0.1, see Section 8 for the evaluation).
Figure 4 demonstrates extremely sparse spiking activity of about $3.3$ Hz in
the recurrent layer (see Fig. 14 for mean firing rate evolution during the
training), nevertheless the network is able to generate the target patterns
well despite PCM non-idealities.

Table 1: Performance evaluation of spiking RNNs with models of PCM crossbar arrays. Method | Sign-gradient | Stochastic | Multi-mem (N=4) | Multi-mem (N=8) | Mixed-precision
---|---|---|---|---|---
MSE Loss | 0.2080 | 0.1808 | 0.1875 | 0.1645 | 0.0380
Figure 4: Dynamics of a network trained with the mixed-precision algorithm.
The raster plot (top) shows the sparse spiking activity of recurrent-layer
neurons. The training loss (bottom left) demonstrates MSE loss over 250 epochs
is averaged over ten best network hyperparameters (see Fig 15 for the best
performing hyperparameter). Properly-tuned neuronal time constants and trained
network weights result in generated patterns following the targets (bottom
right). The generated patterns are extracted from three different spiking
RNNs. Figure 5: (left) The effective conductance distributions $(G^{+}-G^{-})$
of the synapses in the input, recurrent and output layer, at the end of the
training with the mixed-precision method. (right) Averaged over 50 training
runs, the mean number of PCM devices requiring a refresh is shown for each
layer. The refresh operation was not needed for recurrent and output layers.
Figure 6: Total number of WRITE pulses applied to PCM devices are shown for
the input, recurrent and output layers. Only 0.07%, 0.07% and 0.1% of PCM
devices within each layer are programmed respectively during mixed-precision
training.

We also observed that, with our training settings, the weight saturation
problem due to the differential configuration does not occur often and only a
few devices ($<1\%$, as shown in Fig 5 (right)) are required to be refreshed.
This is partially because the total number of WRITE pulses applied during the
training is very low, i.e. only $\sim 12$ WRITE pulses are applied per epoch,
as shown in Fig. 6. Thanks to the mixed-precision algorithm, only large-enough
accumulated gradients lead to the generation of WRITE pulses. Fig. 5 (left)
demonstrates the effective weight distribution of the PCM synapses at the end
of the training.

To simulate an ideal device model, we disabled all noise and drift effects in
the simulation framework (see Supplementary Note 2). We kept a limited weight
resolution of 4-bit, an optimistic but achievable target for PCM devices
Gallo_etal18 . This ideal PCM model is therefore equivalent to a digital 4-bit
memory. Table 2 summarizes the performances of networks utilizing different
memristor-aware weight update methods with this ideal model. The results show
that the stochastic updates, the multi-memristor updates with $N=8$, and the
mixed-precision updates can solve the pattern generation task. Similarly to
the realistic device model results in Table 1, the mixed-precision method
achieved the best accuracy. Unsurprisingly, all methods performed better in
the absence of the previously discussed PCM non-idealities. One remark is that
stochastic updates allow for a better performance than both multi-memristor
methods, indicating the necessity of having a low number of stochastic updates
when the network is trained with the quantized weights.

Moreover, to further evaluate the effect of limited bit precision on the
network performance, we trained the same network with e-prop using standard
FP32 (single-precision floating point) weights. With the the high resolution
of the FP32 training, we achieved an mean squared error (MSE) loss of 0.0215,
which is comparable to mixed-precision training using either ideal quantized
memory or PCM cell model.

Table 2: Performance evaluation of spiking RNNs with an ideal crossbar array model333Multi-memristor configurations are implemented assuming a 4-bit resolution per memory cell. Hence $N=4$ and $N=8$ is equal to having 7-bit and 8-bit resolutions digital weights per synapse, respectively. Method | Sign-gradient | Stochastic | Multi-mem (N=4) | Multi-mem (N=8) | Mixed-precision
---|---|---|---|---|---
MSE Loss | 0.1021 | 0.0758 | 0.1248 | 0.0850 | 0.0289

Similar to Nandakumar_etal18 , we observed that probability scaling factor $p$
in the stochastic update method allows the tuning of how many number of
devices are being programmed during the training. Figure 7 demonstrates that
with the increase of $p$ (decreasing update probability), the number of WRITE
pulses applied to PCM devices can be decreased up to one order of magnitude
without degrading the loss.

Figure 7: The stochastic update method enables tuning the number of WRITE
pulses to be applied to PCM devices.

##  5 Discussion

On-chip training of spiking RNNs enables the low-power deployment of the
intelligent computing systems at the edge with learning and adaptation
capabilities davies_etal21 ; frenkel_etal21 . In this work, we evaluated four
widely used memristor update mechanisms on spiking RNNs based on the ideal
gradient information calculated by the e-prop learning rule. Crossbar
architectures with memristive devices intrinsically support the sparse, event-
driven, and asynchronous processing of spike flow through in-memory computing,
reducing the dynamic power consumption and shrinking the memory area.
Implementing learning on-chip imposes both space and time locality constraints
in recurrent architectures, which can be met by the e-prop learning rule, for
example. By doing extensive hyperparameter exploration, we investigated the
performances of weight update mechanisms despite various PCM non-idealities on
a second-long pattern-generation task. Solving the pattern generation task
demonstrates the recurrent network’s ability to preserve task-relevant
information in its activation dynamics for at least one second period. Such
ability is required for applications processing temporal inputs such as
keyword detection, motor control, and bio-signal processing.

Among the mechanisms we studied, the mixed-precision update leads to the best
accuracy. This is because accumulating instantaneous gradients on high-
precision enables the use of a low learning rate, making the ideal weight
update magnitude comparable to minimum programmable PCM conductance change,
resulting in an improved convergence of the network. Mixed-precision hardware
resembles the concept of the cascade memory model in neuroscience, where
complex synapses hold a hidden state variable which only becomes visible after
hitting a threshold Fusi_Abbott07 . This meta-plastic model has recently been
used for solving catastrophic forgetting in some benchmarks laborieux_etal21 .
However, the mixed-precision scheme comes at the cost of having a high-
precision memory that accumulates the updates. This has previously been done
by introducing a co-processor next to the crossbar memristor array
Nandakumar_etal20a . Yet, accumulating the gradients allows reliable PCM
programming, reducing the number of programming devices and speeds up the
learning. Therefore the synergy between memristor-based synapses and learning
rules that intrinsically allows the accumulation of gradients, e.g., e-prop,
is remarkable and worth exploring.

One of the key levers to increase the learning performance is the memory
resolution, which comes in line with the mixed-precision learning study of
Nandakumar_etal20a . However, this comes at the price of increasing the
synapse size by increasing the number of devices. One solution to this problem
is to employ binary synapses with a stochastic rounding update scheme
Muller_Indiveri15a . Introducing stochasticity reduces the learning rate and,
in expectation, moves the weight parameters to the optimal binary value.
Stochastic rounding approaches have already been employed successfully not
only in memristive devices but also in fully-digital designs to apply low
learning rate values to quantized weights Frenkel_etal19 ; Frenkel_etal20 . It
has also been shown that the intrinsic cycle-to-cycle variability of
memristive devices can implement stochastic updates in an area-efficient
manner Payvand_etal19 .

To the best of our knowledge, there is no report yet on the online training of
spiking RNN with e-prop learning rule based on realistic PCM synapse models.
Our work compares several previously developed methods that are designed to
cope with memristor non-idealities and demonstrates that accumulating
gradients allows more reliable programming of PCM devices, reduces the number
of programming devices and outperforms other synaptic weight-update
mechanisms. Future work will need to evaluate the impacts of the implemented
weight update schemes using more extensive datasets for a more interpretable
benchmarking, further incorporate the PCM devices for emulating temporal
dynamics such as eligibility traces Demirag_etal21a , as well as exploring
other learning rules such as OSTL or Real Time Recurrent Learning (RTRL)
Williams_Zipser89 .

##  6 Methods

For the chosen pattern generation task, the network consists of 100 input and
100 recurrent LIF neurons, and one leaky-integrator output unit. The network
receives a fixed-rate Poisson input, and the target pattern is a one-second-
long sequence defined as the sum of four sinusoids (1 Hz, 2 Hz, 3 Hz and 5
Hz), whose phases and amplitudes are randomly sampled from uniform
distributions $[0,2\pi]$ and $[0.5,2]$, respectively. Throughout the training,
all layer weights {$W_{ji}^{\mathrm{in}}$, $W_{ji}^{\mathrm{rec}}$ and
$W_{kj}^{\mathrm{out}}$}, are kept plastic and the device conductances are
clipped between 0.1 and 12 $\mu$S. This benchmark is adapted from
Bellec_etal18 .

We trained $\sim 1000$ different spiking RNNs for each of the weight update
methods described in Section 2.3. Each network shares the same architecture,
except for their synapse implementations, some of their hyperparameters and
weight initialization. Because each weight update method requires a few
specific additional hyperparameters and may considerably affect the network
dynamics, we tuned the network hyperparameters for each of the update methods
using Bayesian optimization Snoek_etal12a . . Based on network performances
over 250 epochs of the pattern generation task, we selected the best
performing network hyperparameters out of 1000 candidates. By doing so, we
evaluated how well different weight update methods can reflect the ideal
weight update calculated by the e-prop learning rule on a PCM substrate.

##  7 Acknowledgements

This project has received funding from the European Union’s H2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
861153, Memscales project (871371), and ERA-NET CHIST-ERA programme by SNSF
(20CH21186999 / 1).

## References

  * [1] Valentin Khrulkov, Alexander Novikov, and Ivan Oseledets.  Expressive power of recurrent neural networks.  arXiv, 2017.
  * [2] Aaron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu.  WaveNet: A Generative Model for Raw Audio.  arXiv:1609.03499 [cs], 09 2016.
  * [3] Zachary Teed and Jia Deng.  RAFT: Recurrent All-Pairs Field Transforms for Optical Flow.  arXiv, 2020.
  * [4] Tom B Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei.  Language Models are Few-Shot Learners.  arXiv, 2020.
  * [5] Christopher Berner, Greg Brockman, Brooke Chan, Vicki Cheung, Christy Dennison, David Farhi, Quirin Fischer, Shariq Hashme, Chris Hesse, Rafal Józefowicz, Scott Gray, Catherine Olsson, Jakub Pachocki, Michael Petrov, Tim Salimans, Jeremy Schlatter, Jonas Schneider, Szymon Sidor, Ilya Sutskever, Jie Tang, Filip Wolski, and Susan Zhang.  Dota 2 with Large Scale Deep Reinforcement Learning.  page 66, 01 2019.
  * [6] David Ha and Jürgen Schmidhuber.  World Models.  arXiv:1803.10122 [cs, stat], 03 2018.
  * [7] Oriol Vinyals, Igor Babuschkin, Wojciech M. Czarnecki, Michaël Mathieu, Andrew Dudzik, Junyoung Chung, David H. Choi, Richard Powell, Timo Ewalds, Petko Georgiev, Junhyuk Oh, Dan Horgan, Manuel Kroiss, Ivo Danihelka, Aja Huang, Laurent Sifre, Trevor Cai, John P. Agapiou, Max Jaderberg, Alexander S. Vezhnevets, Rémi Leblond, Tobias Pohlen, Valentin Dalibard, David Budden, Yury Sulsky, James Molloy, Tom L. Paine, Caglar Gulcehre, Ziyu Wang, Tobias Pfaff, Yuhuai Wu, Roman Ring, Dani Yogatama, Dario Wünsch, Katrina McKinney, Oliver Smith, Tom Schaul, Timothy Lillicrap, Koray Kavukcuoglu, Demis Hassabis, Chris Apps, and David Silver.  Grandmaster level in StarCraft II using multi-agent reinforcement learning.  Nature, 575(7782):350–354, 10 2019.
  * [8] R.J. Douglas and K.A.C. Martin.  Neocortex.  In G.M. Shepherd, editor, The synaptic organization of the brain, chapter 12, pages 459–509. Oxford University Press, Oxford, New York, 4th edition, 1998.
  * [9] R.J. Douglas, M.A. Mahowald, and C. Mead.  Neuromorphic analogue VLSI.  Annual Review of Neuroscience, 18:255–281, 1995.
  * [10] Friedemann Zenke and Emre O Neftci.  Brain-Inspired Learning on Neuromorphic Substrates.  Proceedings of the IEEE, pages 1–16, 2021.
  * [11] Stefano Ambrogio, Pritish Narayanan, Hsinyu Tsai, Robert M. Shelby, Irem Boybat, Carmelo di Nolfo, Severin Sidler, Massimo Giordano, Martina Bodini, Nathan C. P. Farinha, Benjamin Killeen, Christina Cheng, Yassine Jaoudi, and Geoffrey W. Burr.  Equivalent-accuracy accelerated neural-network training using analogue memory.  Nature, 558(7708):60–67, June 2018.
  * [12] Can Li, Daniel Belkin, Yunning Li, Peng Yan, Miao Hu, Ning Ge, Hao Jiang, Eric Montgomery, Peng Lin, Zhongrui Wang, Wenhao Song, John Paul Strachan, Mark Barnell, Qing Wu, R. Stanley Williams, J. Joshua Yang, and Qiangfei Xia.  Efficient and self-adaptive in-situ learning in multilayer memristor neural network.  Nature Communications, 9(2385):1–8, June 2018.
  * [13] Mirko Prezioso, Farnood Merrikh-Bayat, BD Hoskins, GC Adam, Konstantin K Likharev, and Dmitri B Strukov.  Training and operation of an integrated neuromorphic network based on metal-oxide memristors.  Nature, 521(7550):61–64, 2015.
  * [14] Thomas Dalgaty, Niccolo Castellani, Clément Turck, Kamel-Eddine Harabi, Damien Querlioz, and Elisa Vianello.  In situ learning using intrinsic memristor variability via markov chain monte carlo sampling.  Nature Electronics, 4(2):151–161, 2021.
  * [15] Fuxi Cai, Suhas Kumar, Thomas Van Vaerenbergh, Xia Sheng, Rui Liu, Can Li, Zhan Liu, Martin Foltin, Shimeng Yu, Qiangfei Xia, et al.  Power-efficient combinatorial optimization using intrinsic noise in memristor hopfield neural networks.  Nature Electronics, 3(7):409–418, 2020.
  * [16] Abu Sebastian, Manuel Le Gallo, and Evangelos Eleftheriou.  Computational phase-change memory: beyond von Neumann computing.  Journal of Physics D: Applied Physics, 52(44):443002, August 2019\.
  * [17] J. Backus.  Can programming be liberated from the von Neumann style?: a functional style and its algebra of programs.  Communications of the ACM, 21(8):613–641, 1978.
  * [18] G. Indiveri and S.-C. Liu.  Memory and information processing in neuromorphic systems.  Proceedings of the IEEE, 103(8):1379–1397, 2015.
  * [19] Melika Payvand, Manu V Nair, Lorenz K Müller, and Giacomo Indiveri.  A neuromorphic systems approach to in-memory computing with non-ideal memristive devices: From mitigation to exploitation.  Faraday Discussions, 213:487–510, 2019.
  * [20] Elisabetta Chicca and Giacomo Indiveri.  A recipe for creating ideal hybrid memristive-CMOS neuromorphic processing systems.  Applied Physics Letters, 116(12):120501, 2020.
  * [21] Xiaochen Peng, Shanshi Huang, Yandong Luo, Xiaoyu Sun, and Shimeng Yu.  DNN+NeuroSim: An End-to-End Benchmarking Framework for Compute -in-Memory Accelerators with Versatile Device Technologies.  IEEE International Electron Devices Meeting (IEDM), pages 32.5.1–32.5.4, 2019.
  * [22] Xiaochen Peng, Shanshi Huang, Hongwu Jiang, Anni Lu, and Shimeng Yu.  DNN+NeuroSim V2.0: An End-to-End Benchmarking Framework for Compute-in-Memory Accelerators for On-chip Training.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, PP(99):1–1, 2020.
  * [23] Geoffrey W. Burr, Matthew J. Brightsky, Abu Sebastian, Huai-Yu Cheng, Jau-Yi Wu, Sangbum Kim, Norma E. Sosa, Nikolaos Papandreou, Hsiang-Lan Lung, Haralampos Pozidis, Evangelos Eleftheriou, and Chung H. Lam.  Recent Progress in Phase-Change Memory Technology.  IEEE Journal on Emerging and Selected Topics in Circuits and Systems, 6(2):146–162, 06 2016.
  * [24] Geoffrey W Burr, Robert M Shelby, Abu Sebastian, Sangbum Kim, Seyoung Kim, Severin Sidler, Kumar Virwani, Masatoshi Ishii, Pritish Narayanan, Alessandro Fumarola, Lucas L Sanches, Irem Boybat, Manuel Le Gallo, Kibong Moon, Jiyoo Woo, Hyunsang Hwang, and Yusuf Leblebici.  Neuromorphic computing using non-volatile memory.  Advances in Physics: X, 2(1):89–124, 2017.
  * [25] Tomas Tuma, Angeliki Pantazi, Manuel Le Gallo, Abu Sebastian, and Evangelos Eleftheriou.  Stochastic phase-change neurons.  Nature Nanotechnology, 11(8):693–699, 2016.
  * [26] Geethan Karunaratne, Manuel Le Gallo, Giovanni Cherubini, Luca Benini, Abbas Rahimi, and Abu Sebastian.  In-memory hyperdimensional computing.  Nature Electronics, 3(6):327–337, 06 2020.
  * [27] Yigit Demirag, Filippo Moro, Thomas Dalgaty, Gabriele Navarro, Charlotte Frenkel, Giacomo Indiveri, Elisa Vianello, and Melika Payvand.  PCM-trace: Scalable Synaptic Eligibility Traces with Resistivity Drift of Phase-Change Materials.  2021 IEEE International Symposium on Circuits and Systems (ISCAS), pages 1–5, 2021.
  * [28] Manuel Le Gallo, Abu Sebastian, Giovanni Cherubini, Heiner Giefers, and Evangelos Eleftheriou.  Compressed Sensing With Approximate Message Passing Using In-Memory Computing.  IEEE Transactions on Electron Devices, 65(10):4304–4312, 2018.
  * [29] SR Nandakumar, Manuel Le Gallo, Irem Boybat, Bipin Rajendran, Abu Sebastian, and Evangelos Eleftheriou.  A phase-change memory model for neuromorphic computing.  Journal of Applied Physics, 124(15):152135, 2018.
  * [30] Guillaume Bellec, Franz Scherr, Anand Subramoney, Elias Hajek, Darjan Salaj, Robert Legenstein, and Wolfgang Maass.  A solution to the learning dilemma for recurrent networks of spiking neurons.  Nature Communications, 11(3625):1–15, 2020.
  * [31] Yigit Demirag.  Multiphysics modeling of $\text{Ge}_{2}\text{Sb}_{2}\text{Te}_{5}$ based synaptic devices for brain inspired computing.  Master’s thesis, Ihsan Dogramaci Bilkent University, Ankara, Turkey, July 2018.
  * [32] Manuel Le Gallo and Abu Sebastian.  An overview of phase-change memory device physics.  Journal of Physics D: Applied Physics, 53(21):213002, 2020.
  * [33] Manuel Le Gallo, Matthias Kaes, Abu Sebastian, and Daniel Krebs.  Subthreshold electrical transport in amorphous phase-change materials.  New Journal of Physics, 17(9):093035, 09 2015.
  * [34] Manuel Le Gallo, Aravinthan Athmanathan, Daniel Krebs, and Abu Sebastian.  Evidence for thermally assisted threshold switching behavior in nanoscale phase-change memory cells.  Journal of Applied Physics, 119(2):025704, 01 2016.
  * [35] D. Ielmini, S. Lavizzari, D. Sharma, and A. L. Lacaita.  Physical Interpretation, Modeling and Impact on Phase Change Memory (PCM) Reliability of Resistance Drift Due to Chalcogenide Structural Relaxation.  2007 IEEE International Electron Devices Meeting, pages 939–942, 2007.
  * [36] IV Karpov, M Mitra, D Kau, G Spadini, YA Kryukov, and VG Karpov.  Fundamental drift of parameters in chalcogenide phase change memory.  Journal of Applied Physics, 102(12):124503, 2007.
  * [37] A. Redaelli, A. Pirovano, A. Benvenuti, and A. L. Lacaita.  Threshold switching and phase transition numerical models for phase change memory simulations.  Journal of Applied Physics, 103(11):111101, 2008.
  * [38] Martin Salinga, Egidio Carria, Andreas Kaldenbach, Manuel Bornhöfft, Julia Benke, Joachim Mayer, and Matthias Wuttig.  Measurement of crystal growth velocity in a melt-quenched phase-change material.  Nature Communications, 4(1):2371, 2013.
  * [39] M Nardone, V I Kozub, I V Karpov, and V G Karpov.  Possible mechanisms for 1/f noise in chalcogenide glasses: A theoretical description.  Physical Review B, 79(16):165206, 2009.
  * [40] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [41] Blake A. Richards, Timothy P. Lillicrap, Philippe Beaudoin, Yoshua Bengio, Rafal Bogacz, Amelia Christensen, Claudia Clopath, Rui Ponte Costa, Archy de Berker, Surya Ganguli, Colleen J. Gillon, Danijar Hafner, Adam Kepecs, Nikolaus Kriegeskorte, Peter Latham, Grace W. Lindsay, Kenneth D. Miller, Richard Naud, Christopher C. Pack, Panayiota Poirazi, Pieter Roelfsema, João Sacramento, Andrew Saxe, Benjamin Scellier, Anna C. Schapiro, Walter Senn, Greg Wayne, Daniel Yamins, Friedemann Zenke, Joel Zylberberg, Denis Therien, and Konrad P. Kording.  A deep learning framework for neuroscience.  Nature Neuroscience, 22(11):1761–1770, 11 2019.
  * [42] Nicolas Frémaux and Wulfram Gerstner.  Neuromodulated spike-timing-dependent plasticity, and theory of three-factor learning rules.  Front. Neur. Circ., 9:85, 2016.
  * [43] João Sacramento, Rui Ponte Costa, Yoshua Bengio, and Walter Senn.  Dendritic cortical microcircuits approximate the backpropagation algorithm.  In Advances in Neural Information Processing Systems, pages 8721–8732, 2018.
  * [44] Isabella Pozzi, Sander Bohté, and Pieter Roelfsema.  A Biologically Plausible Learning Rule for Deep Learning in the Brain.  arXiv, 2018.
  * [45] Wulfram Gerstner and Werner M Kistler.  Spiking neuron models: Single neurons, populations, plasticity.  Cambridge university press, Cambridge, United Kingdom, 2002.
  * [46] P. Dayan and L.F. Abbott.  Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems.  MIT Press, Cambridge, MA, USA, 2001.
  * [47] D. Sussillo and L.F. Abbott.  Generating coherent patterns of activity from chaotic neural networks.  Neuron, 63(4):544–557, 2009.
  * [48] Peter U Diehl, Daniel Neil, Jonathan Binas, Matthew Cook, Shih-Chii Liu, and Michael Pfeiffer.  Fast-classifying, high-accuracy spiking deep networks through weight and threshold balancing.  In International Joint Conference on Neural Networks (IJCNN), pages 1–8. IEEE, 2015.
  * [49] Wilten Nicola and Claudia Clopath.  Supervised Learning in Spiking Neural Networks with FORCE Training.  Nature Communications, 8(1):2208, 12 2017.
  * [50] Emre O Neftci, Hesham Mostafa, and Friedemann Zenke.  Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks.  IEEE Signal Processing Magazine, 36(6):51–63, 2019.
  * [51] Jun Haeng Lee, Tobi Delbruck, and Michael Pfeiffer.  Training deep spiking neural networks using backpropagation.  Frontiers in Neuroscience, 10:508, 2016.
  * [52] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [53] Thomas Bohnstingl, Stanisław Woźniak, Wolfgang Maass, Angeliki Pantazi, and Evangelos Eleftheriou.  Online Spatio-Temporal Learning in Deep Neural Networks.  arXiv, 2020.
  * [54] Ning Qiao, Hesham Mostafa, Federico Corradi, Marc Osswald, Fabio Stefanini, Dora Sumislawska, and Giacomo Indiveri.  A reconfigurable on-line learning spiking neuromorphic processor comprising 256 neurons and 128k synapses.  Frontiers in neuroscience, 9:141, 2015.
  * [55] Melika Payvand, Lorenz K Muller, and Giacomo Indiveri.  Event-based circuits for controlling stochastic learning with memristive devices in neuromorphic architectures.  In Circuits and Systems (ISCAS), 2018 IEEE International Symposium on, pages 1–5. IEEE, 2018.
  * [56] Melika Payvand, Mohammed E. Fouda, Fadi Kurdahi, Ahmed M. Eltawil, and Emre O. Neftci.  On-Chip Error-Triggered Learning of Multi-Layer Memristive Spiking Neural Networks.  IEEE Journal on Emerging and Selected Topics in Circuits and Systems, 10(4):522–535, 2020.
  * [57] Irem Boybat, Manuel Le Gallo, Timoleon Moraitis, Thomas Parnell, Tomas Tuma, Bipin Rajendran, Yusuf Leblebici, Abu Sebastian, and Evangelos Eleftheriou.  Neuromorphic computing with multi-memristive synapses.  Nature Communications, 9:2514, 2018.
  * [58] S. R. Nandakumar, Manuel Le Gallo, Christophe Piveteau, Vinay Joshi, Giovanni Mariani, Irem Boybat, Geethan Karunaratne, Riduan Khaddam-Aljameh, Urs Egger, Anastasios Petropoulos, Theodore Antonakopoulos, Bipin Rajendran, Abu Sebastian, and Evangelos Eleftheriou.  Mixed-Precision Deep Learning Based on Computational Memory.  Frontiers in Neuroscience, 14:406, 2020.
  * [59] Manu V Nair, Lorenz K. Mueller, and Giacomo Indiveri.  A differential memristive synapse circuit for on-line learning in neuromorphic computing systems.  Nano Futures, 1(3):1–12, 2017.
  * [60] G.W. Burr, R.M. Shelby, C.di Nolfo, J.W. Jang, R.S. Shenoy, P. Narayanan, K. Virwani, E.U. Giacometti, B. Kurdi, and H. Hwang.  Experimental Demonstration and Tolerancing of a Large-Scale Neural Network (165,000 Synapses), using Phase-Change Memory as the Synaptic Weight Element.  2014 IEEE International Electron Devices Meeting, pages 1–4, 2014\.
  * [61] Lukas Balles, Fabian Pedregosa, and Nicolas Le Roux.  The Geometry of Sign Gradient Descent.  arXiv, 2020.
  * [62] Manu V Nair and Piotr Dudek.  Gradient-descent-based learning in memristive crossbar arrays.  In International Joint Conference on Neural Networks (IJCNN), pages 1–7. IEEE, 2015.
  * [63] L. Müller, M. Nair, and G. Indiveri.  Randomized unregulated step descent for limited precision synaptic elements.  In International Symposium on Circuits and Systems, (ISCAS). IEEE, 2017.
  * [64] Melika Payvand, Mohammed E. Fouda, Fadi Kurdahi, Ahmed M. Eltawil, and Emre O. Neftci.  On-Chip Error-Triggered Learning of Multi-Layer Memristive Spiking Neural Networks.  IEEE Journal on Emerging and Selected Topics in Circuits and Systems, 10(4):522–535, 2020.
  * [65] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [66] Aravinthan Athmanathan, Milos Stanisavljevic, Nikolaos Papandreou, Haralampos Pozidis, and Evangelos Eleftheriou.  Multilevel-Cell Phase-Change Memory: A Viable Technology.  IEEE Journal on Emerging and Selected Topics in Circuits and Systems, 6(1):87–100, 2016.
  * [67] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio.  Quantized Neural Networks: Training Neural Networks with Low Precision Weights and Activations.  arXiv, 2016.
  * [68] Mike Davies, Andreas Wild, Garrick Orchard, Yulia Sandamirskaya, Gabriel A Fonseca Guerra, Prasad Joshi, Philipp Plank, and Sumedh R Risbud.  Advancing neuromorphic computing with loihi: A survey of results and outlook.  Proceedings of the IEEE, 109(5):911–934, 2021.
  * [69] Charlotte Frenkel, David Bol, and Giacomo Indiveri.  Bottom-up and top-down neural processing systems design: Neuromorphic intelligence as the convergence of natural and artificial intelligence.  arXiv preprint arXiv:2106.01288, 2021.
  * [70] S. Fusi and L.F. Abbott.  Limits on the memory storage capacity of bounded synapses.  Nature Neuroscience, 10:485–493, 2007.
  * [71] Axel Laborieux, Maxence Ernoult, Tifenn Hirtzlin, and Damien Querlioz.  Synaptic metaplasticity in binarized neural networks.  Nature communications, 12(1):1–12, 2021.
  * [72] Lorenz K Muller and Giacomo Indiveri.  Rounding methods for neural networks with low resolution synaptic weights.  arXiv preprint arXiv:1504.05767, pages 1–11, 2015.
  * [73] Charlotte Frenkel, Jean-Didier Legat, and David Bol.  Morphic: A 65-nm 738k-synapse/mm2 quad-core binary-weight digital neuromorphic processor with stochastic spike-driven online learning.  IEEE Transactions on Biomedical Circuits and Systems, 13(5):999–1010, 2019.
  * [74] Charlotte Frenkel, Jean-Didier Legat, and David Bol.  A 28-nm convolutional neuromorphic processor enabling online learning with spike-based retinas.  In 2020 IEEE International Symposium on Circuits and Systems (ISCAS), pages 1–5. IEEE, 2020.
  * [75] Ronald J. Williams and David Zipser.  A Learning Algorithm for Continually Running Fully Recurrent Neural Networks.  Neural Computation, 1(2):270–280, 1989.
  * [76] Guillaume Bellec, Darjan Salaj, Anand Subramoney, Robert Legenstein, and Wolfgang Maass.  Long short-term memory and learning-to-learn in networks of spiking neurons.  In Advances in Neural Information Processing Systems, pages 787–797, 2018.
  * [77] Jasper Snoek, Hugo Larochelle, and Ryan P. Adams.  Practical bayesian optimization of machine learning algorithms.  In Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 2, NIPS’12, page 2951–2959, Red Hook, NY, USA, 2012. Curran Associates Inc.

##  8 Supplementary

### Supplementary Note 1

We implemented the PCM crossbar array simulation framework in PyTorch [40],
which can be used for both the inference and the training of Artificial Neural
Networks or SNNs. Built on top of the statistical model introduced by
Nandakumar et al. [29], our crossbar model supports asynchronous SET, RESET
and READ operations over entire crossbar structures and simultaneously keep
tracks of the temporal evolution of device conductances.

A single crossbar array consists of $P\times Q$ nodes (each node representing
a synapse), where single node has $2N$ memristors arranged using the
differential architecture ($N$ potentiation, $N$ depression devices). Each
memristor state is represented by four variables, $t_{p}$ for storing the last
time the device is written (which is used to calculate the effect of the
drift), $count$ for counting how many times it has written (to be used later
in the arbiter of N-memristor architectures), $P_{mem}$ for its programming
history (required by the PCM model) and $G$ for representing the conductance
of the the device at $T_{0}$ seconds later after the last programming time.
The initial conductances of PCM devices in the crossbar array are assumed to
be iterativelly programmed to High-Resistive State (HRS), sampled from a
Normal distribution $\mathcal{N}(\mu=0.1,\sigma=0.01)$ $\mu$S.

The PCM crossbar simulation framework supports three major functions: READ,
SET and RESET. The READ function takes the pulse time of the applied READ
pulse, $t$, and calculates the effect of drift based on the last programming
time $t_{p}$. Then, it adds the conductance-dependent READ noise and returns
the conductance values of whole array. The SET function takes the timing
information of the applied SET pulse, together with a mask of shape $(2\times
N\times P\times Q)$ and calculates the effect of the application of a single
SET pulse on the PCM devices that are selected with the mask. Finally, the
RESET function initializes all the state variables of devices selected with a
mask and initializes the conductances using a Normal distribution
$\mathcal{N}(\mu=0.1,\sigma=0.01)$ $\mu$S.

### Supplementary Note 2

READ and WRITE operations to simulated PCM devices in the crossbar model are
stochastic and subject to the temporal conductance drift. Additionally, PCM
devices offer a very limited bit precision. Therefore, to ease the network
training procedure, especially the hyperparameter tuning, we developed the
perf-mode. When crossbar model is operated in perf-mode, all stochasticity
sources and the conductance drift are disabled. READ operations directly
access the device conductance without $1/f$ noise and drift, whereas SET
operations increase the device conductance as

$G_{N}=G_{N-1}+\frac{G_{MAX}}{2^{CB_{RES}}},$ (4)

where, $G_{MAX}$ is the maximum PCM conductance set to 12 $\mu$S (conductivity
boundaries are determined based on the device measurements from [29]), and
$CB_{RES}$ is the desired bit-resolution of a single PCM device. In a
nutshell, the perf-mode turns PCMs into an ideal memory cells corresponding to
a digital memory with a limited bit precision.

(a) Comparison of the full PCM model and its
perf-mode equivalent after 8 consecutive SET
pulses.

(b) Comparison of the full PCM model and its
perf-mode equivalent after 8 consecutive SET
pulses, averaged over 300 measurements showing
the effect of drift.

Figure 8: The PCM crossbar model supports both the full PCM model from [29]
and its corresponding simplified version as an ideal digital memory in perf-
mode.

### Supplementary Note 3

Here, we demonstrate the impact of using multiple memristor devices per
synapse (arranged in differential configuration) on the precision of targeted
programming updates. Specifically, we modeled synapses with $N={1,4,8}$ PCM
devices and programmed them from initial conditions of integer conductance
values $G_{source}\in\\{-10,10\\}$ $\mu$ S to integer conductance values
$G_{target}\in\\{-10,10\\}$ $\mu$ S using the multi-memristor update scheme
described in Section 2.3. The effective conductance of a synapse is calculated
by $G_{syn}=\sum_{i=0}^{N}G_{i}^{+}-\sum_{i=0}^{N}G_{i}^{-}$, however we
normalized the conductance across 1-PCM, 4-PCM and 8-PCM architectures for an
easier comparison, such that
$G_{syn}=\frac{1}{N}(\sum_{i=0}^{N}G_{i}^{+}-\sum_{i=0}^{N}G_{i}^{-})$.

Our empirical results verifies the claim of Boybat et al. [57] that the
standard deviation and the update resolution of the write process decreases by
$\sqrt{N}$.

Figure 9: Multi-memristor configuration with 1 PCM (one depression and one
potentiation) per synapse Figure 10: Multi-memristor configuration with 8 PCM
(four depression and four potentiation) per synapse Figure 11: Multi-memristor
configuration with 16 PCM (eight depression and eight potentiation) per
synapse

### Supplementary Note 4

In the differential architectures, consecutive SET pulses applied to positive
and negative memristors may cause the saturation of the synaptic conductance
and block further updates. The saturation effect is more apparent when a
single synapse gets 10+ updates in one direction (potentiation or depression)
during the training. For example, this effect is clearly visible in Fig. 9,
Fig. 10 and Fig. 11, when the source conductance and target conductances
differ by more than 8-10 $\mu$S.

We implemented a weight update scheme denoted as the update-ready criterion,
which aims to prevent conductance saturation while applying single large
updates. Before doing the update, we read both positive and negative pair
conductances, and check if the target update is possible. If not, we reset
both devices, calculate the new target and apply the number of pulses
accordingly. For example, given $G^{+}=8\mu$S and $G^{-}=4\mu$S and the
targeted update $+6\mu$S, the algorithm decides to reset both devices because
$G^{+}$ can’t be increased up to $14\mu$S. After both devices are reset,
$G^{+}$ can be programmed to $10\mu$S). Although our PCM crossbar array
simulation framework supports it, this weight transfer criterion is not used
in our simulations because it requires reading the device states during the
update.

Figure 12: Update-ready criterion tested with $N=1$ memristor per synapse.

### Supplementary Note 5

We have defined the task success criteria as MSE Loss $<0.1$ based on visual
inspection. Below in Fig 13, some network performances are shown.

Figure 13: Comparison of network performances with six different loss values.

### Supplementary Note 6

Figure 14: Mean firing rate of 50 networks with PCM synapses trained using the
mixed-precision method. Figure 15: MSE loss of 50 networks trained with PCM
synapses using the mixed-precision method.

  *[RNNs]: Recurrent Neural Network
  *[PCM]: Phase Change Memory
  *[RNN]: Recurrent Neural Network
  *[RL]: Reinforcement Learning
  *[SNN]: Spiking Neural Network
  *[SNNs]: Spiking Neural Network
  *[CMOS]: Complementary Metal-Oxide-Semiconductor
  *[SRAM]: Static Random Access Memory
  *[PCMs]: Phase Change Memory
  *[GST]: Ge2Sb2Te5
  *[LIF]: Leaky Integrate-and-Fire
  *[BPTT]: Back-propagation Through Time
  *[LSTM]: Long Short Term Memory
  *[OSTL]: Online Spatio-Temporal Learning
  *[LB]: Learning Block
  *[RTRL]: Real Time Recurrent Learning
  *[HRS]: High-Resistive State