# The Challenges Ahead for Bio-inspired Neuromorphic Event Processors: How
Memristors Dynamic Properties Could Revolutionize Machine Learning

Marco Rasetto
University of Pittsburgh, Center for Neural Basis of Cognition
Pittsburgh, USA
<EMAIL_ADDRESS>
&Qingzhou Wan
University of Pittsburgh,
Pittsburgh, USA
&Himanshu Akolkar
University of Pittsburgh, Center for Neural Basis of Cognition
Pittsburgh, USA
&Bertram Shi
Hong Kong University of Science and Technology,
Hong Kong, China
&Feng Xiong
University of Pittsburgh,
Pittsburgh, USA
&Ryad Benosman
University of Pittsburgh, Carnegie Mellon Institute
Pittsburgh, USA
<EMAIL_ADDRESS>
###### Abstract

Neuromorphic engineering has led to the necessary process of rethinking of how
we process and integrate information, analyze data, and use the resulting
insights to improve computation and avoid the current high power and latency
of Artificial Intelligence (AI) hardware. Current neuromorphic processors are,
however, limited by digital technologies, which cannot reproduce the abilities
of biological neural computation in terms of power, latency and area cost. In
this paper, we show that the combined use of the dynamic properties of
memristors to implement a model of synaptic integration and the determination
of the correct level of abstraction of biological neural networks has the
potential to open a new range of capabilities for neuromorphic processors. We
test this approach using a novel three-terminal
Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ electrochemical memristor, by deriving
its conductance model and using it to emulate synaptic temporal kernel
computation in the context of a pattern recognition task. We show that these
devices allow for robust results with no loss in precision while opening the
path for an energy efficient approach to build novel bio-inspired processing
units in silicon.

##  1 Introduction

The last decade brought considerable progress in AI, allowed by the advent of
GPUs and similar hardware. However, this progress has not been matched from
the perspective of emulating general intelligence and cognition. Ideas such as
deep multilayer learning has helped solve a certain class of very well defined
problems, but require high energy and large amounts of labeled data [1]. In
addition, important functionalities, such as continuous always-on learning
with a reasonable power budget, are still out of reach. This limits on-board
”intelligence” and reduces autonomy drastically.

Neuromorphic engineering holds the promise of mitigating these restrictions
[2, 3]. Recently, this field has reached a level of maturity that allows it to
impact several other domains where autonomy and low power on-the-edge
computation are crucial. Event cameras are among the most visible development
of neuromorphic engineering [4, 5]. These sensors are revolutionizing
perception, especially for scenarios involving highly unpredictable dynamic
environments [6]. Although they are perfect candidates to solve a variety of
dynamic robotics perception tasks, their use is currently slowed down by the
absence of an electronic brain-inspired processor, able to process their
output in general computational scenarios.

Currently available neuromorphic processors are still in their infancy. They
are built on the concept of replicating biological neurons into silicon [7, 8,
9, 10]. For number of reasons, their use cases have been limited to specific
narrow scenarios. The most obvious reason is that although we know much about
the brain, we are still missing a complete theory explaining how biological
brains operate. Another reason is that silicon and biological brains are
implemented on different physical substrates, with different capabilities and
constraints.

The very original definition of neuromorphic computing can be summarized as
developing analog computing hardware organized similarly to the brain and that
is able to learn. Today, three decades after that initial movement, it appears
legitimate to wonder how closely should neuromorphic computing hardware be
targeted to emulate brain functionalities? And more importantly, how much of
the biological properties should we aim to replicate in silicon? The baseline
of abstraction being the relation between current back-propagation based
artificial neural networks and biological learning. It seems essential that
the scale of success of this endeavour relies on the ability to achieve
extreme energy efficiency while allowing novel cognitive processes able to
compete and outperform against brute-force conventional GPU back-propagation
based AI.

We argue in this paper that the solution should not rely on a systematic
mimicking of brain features, the reason being that millions of years of
evolution constrained by the biological brain biomaterials ”susbtrate” is the
reason of the current architecture of brains.

###  1.1 Hybrid Substrates: Modeling the correct abstraction level

Effectively exploring biological inspiration does not imply directly copying
nature. Rather, the goal should be to understand the core principles and
abstract them in a way that makes them amenable to implement in available
computational substrates. Today’s technology allows a much wider range of
available materials that could be a better match to design improved silicon
brains based on a different and perhaps more advanced set of computation
building blocks and design rules. Science is full of examples showing the
importance of understanding compared to replicating. Perhaps the most known
one relates to the early days of aviation. Although birds flap wings to fly,
available materials at the beginning of the last century such as wood and
canvas did not have the structural properties to sustain bird-like flight.
These, however, proved to be adequate in designing fixed wing airplanes thanks
to an understanding of the laws of aerodynamics.

The same principle applies to brain-inspired computing. The current open
question is to determine the adequate level of abstraction [11]. There is a
close coupling between the level of abstraction, the physical substrate and
the computational model design. The commonly used level of abstraction, that
is perhaps closest to copying biology directly, models neural computation
using a set of coupled ordinary differential equations. The temporal dynamics
enable the model to integrate information over time. The coupling across state
variables models the integration of information over space [12]. This level of
abstraction is most commonly implemented using analog continuous time VLSI
circuits [13]. These typically offer low power, but suffer from drawbacks
including mismatch and limited programmability.

A closely related level of abstraction uses a set of coupled difference
equations. Models at this level of abstraction are commonly implemented using
digital design methodologies in standard CMOS processes or FPGAs [7, 8, 9,
10]. While they can be implemented using standard CPUs, the timescales of
operation are slower than those of biological neurons. Updates occur at the
level of changes in the input, e.g. frame rates of 30 fps. Examples of models
implemented in this way include most deep learning models, e.g. CNNs [14],
LSTMs [15], RNNs [16], etc. However, even with fully custom designs, these do
not reach the low levels of power sought by neuromorphic engineers.

Hybrid substrates offer an intriguing design space, which has the potential to
leverage the unique advantages of the analog and digital worlds [17]. However,
the level of abstraction best suited for a hybrid implementation is not yet
clear. The most common hybrid model uses analog circuits for computation and
digital circuits for communication [18, 19, 20], leveraging the insight that
the speed of digital circuits is much higher than the typical spike rate of
neurons. Thus, a single digital bus can be multiplexed to carry signals from
multiple neurons. The dominant implementation of this level of abstraction are
multi-chip Address Event Representation (AER) networks [21, 22]. Computation
within each chip is done using analog continuous time circuits. Communication
among these circuits is digital and often asynchronous.

Here we argue in favor of such a hybrid model but using an entirely different
partitioning in the abstraction. Rather than splitting along function
(computation vs communication), we propose a split along dimension (time
versus space). We argue that temporal integration of signals, especially
spiking signals, is best implemented in analog, whereas spatial integration is
best implemented in digital. In particular, newly developed memristive
technologies [23, 24, 25, 26] offer an excellent physical substrate for
temporal integration. In contrast to temporal integration, spatial
integration, which requires communication of signals across space, is best
suited for digital technologies.

###  1.2 Memristors as temporal computation units

Our proposal runs counter to the most dominant paradigm for neuromorphic
memristive systems, which follows the computation/communication partitioning.
Typically, proposals for memristive neuromorphic systems use the memristors in
the communication layer, due to the ability to program their resistance via
current pulses[27, 28, 29, 30]. For example, the commonly proposed memristive
crossbar array uses memristors to implement both weight storage and weight
combination. Weights are stored as the conductances of the memristors, which
are connected so that currents can be summed at a common node.

Following this paradigm, memristors over the decades have found use in many
network-based computational approaches. Doygu et al. simulated memristor
networks able to learn sequences of inputs [31]. Suri et al. used another
simulated network of phase-changing memristors to learn MNIST letters[32].
Systems that perform STDP (Spike Time Dependent Plasticity) using RRAM
(REsistive RAM) devices [33] or that implement recursive networks using PCM
memristors have also been reported [34]. Vicent et al. simulated a network of
STT-MRAM devices for car detection [35]. Networks of memristors have also been
proposed to implement the k-means algorithm [36] and unsupervised learning
[37]. In all of these examples, the synaptic devices modeled using memristors
work in a “static” fashion, i.e. as a fixed scalar multiplier of spike events,
before integration by the “neuron membrane”. This approach has been widely
adopted because of its simplicity of implementation, and its mathematical
tractability, as synaptic operations can be described by simple linear
algebra.

We argue here that volatile memristors exhibiting spontaneous decay,
reminiscent of the decay of excitatory post-synaptic potentials (EPSPs)[38],
are better suited to serve as excellent substrates for dynamic temporal
integration, especially of signals encoded as spikes. This is in part because
these decays are an intrinsic outcome of the device physics and not a circuit-
level design function. Thus, memristors offer an extremely compact (single
device) and energy efficient implementation of the most critical operation in
dynamic computation circuits: the temporal integration [39].

Memristors are ideal in this respect because in real biological networks,
synapses exhibit a far more nuanced repertoire of complex behaviors, with
temporal dynamics, adaptation and short term effects [40]. For instance, when
responding to a single pulse, chemical synapses release neurotransmitters that
can be modelled using exponential decay functions with time constants in the
order of milliseconds or longer [38]. When combined in multiple dendritic
compartments, these effects can encode complex spatio-temporal patterns that
can be learned and recognized by single neurons [41, 42]. Moreover, short term
plasticity (STP) allows synapses to act as gain control, adapt stimuli
selectivity, select temporal frequencies [43], improve temporal stimuli
recognition [44, 45] and enrich neural dynamics [46, 47]. Interestingly, both
STP and exponential decay dynamics have been reproduced in different types of
memristive devices in the past, from atomic switch CBRAM [23], to filamentary
RRAM [24], on electrochemical [25] and on 2D material devices [26]. Memristors
make it feasible to implement more complex synaptic models using single
devices. Although a few papers have explored the computational advantages of
the dynamic properties of memristors, none have proposed a full architecture
capable of fully leveraging the capabilities of memristor dynamics and tested
it on a real-world problem, as described here. Rather, they have focused on
small components. For example, Chang et al. [48] and Lim et al. [49] proposed
that STP could be used to implement a form of short term memory. Berdan et al.
[50] and Wan et al. [45] demonstrated their potential use for limited forms of
spatio-temporal and temporal pattern recognition.

###  1.3 A novel memristor for spatio-temporal computation

Our proposed partitioning is reminiscent of space-time separability:
integration of information over time first, pixel by pixel or more generally,
neuron by neuron, followed by spatial integration of signals across
pixels/neurons. It is well known that space-time separability in digital
signal processing algorithms, offers compelling implementation advantages
[51].

Unfortunately, not all video processing algorithms are space time separable.
The same is true of neuromorphic algorithms. For example, it is certainly not
true that any set of coupled differential equations can be expressed as a
space-time separable set of operations. Indeed, the vast majority cannot.

However, we demonstrate here that a very large and compelling class of
algorithms able to operate in an abstract scale not requiring spiking neurons,
namely, those using time surfaces and hierarchies of time surfaces (HOTS)[52,
53, 54] are by their very nature, space-time separable in the manner described
above, and are ideally suited for implementation using a combination of
memristive devices for temporal integration followed by digital spatial
integration (in particular clustering and mapping to nearest cluster centers).
We combine HOTS and a novel three terminal memristor -
Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ to demonstrate that such architectures
can be used for pattern recognition while being robust to the types of
nonidealities (i.e. random mismatch and noise) encountered in memristive
devices.

###  1.4 Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ memristor - a three terminal
memristor

In recent years, the development of new memristive technologies has allowed
temporal computing and modeling of temporal kernels and short term plasticity
to be carried out on a smaller footprint without the requirements of
additional circuitry [23, 24, 25, 26]. This feature of memristors with
“volatile” properties is extremely promising for neuromorphic networks that
need physical substrate capable of temporal computation that could easily be
scaled to adapt for multiple architectures and algorithms [50]. Among these
devices, electrochemical memristors [55] have become good candidates as
“synaptic devices” thanks to low power, linear and symmetric response, low
variations, and high reliability [55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65].
We have previously proposed the use of a novel electrochemical memristor [45]
based on Lithium Ions and Tungsten Oxide
(Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$) where the conductance can be
programmed as exponential decays in the order of tens to hundred milliseconds
which are optimal for temporal integration of event based stimuli [52, 38].
This device (Fig. 1a) is built by deposition of WO${}_{\textbf{3}}$ films
acting as gate and channel (red) on a LaALO${}_{\textbf{3}}$ substrate (green)
and connected to gold terminals (yellow), forming a three terminal structure.
Thanks to its structure, the Tungsten Oxide can embed additional Lithium ions
that act as short or long term doping charge by changing conductance of the
material [45]. The lithium ions are introduced to the substrate by the
application of a drop of electrolyte gel (blue) on top of the device. In the
next section, we will present a deeper insight on the dynamics of operation of
the Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ memristor and present an improved
computational model that allow us to simulate the device behaviour and
stochasticity, and test it within a full neuromorphic architecture.

##  2 Modeling the electrochemical Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$
memristor dynamics

Unlike the conventional two-terminal memristors, the
Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ electrochemical synapse [Fig. 1] or
electrochemical random access memory (ECRAM) [66], is composed of three
terminals [(S)ource, (D)rain, (G)ate in Fig. 1(a)], similar to a transistor
structure. The conductance between the source and the drain terminal can be
considered as the synaptic weight of the device, which can be read at any
moment by applying a small DC reading voltage (VDS = 0.1 V) across the source
and drain. The resulting current across the two terminals is then measured to
calculate the conductance, GDS, given by GDS = IDS/VDS. This conductance can
be modulated or by applying a square voltage “write” pulse between the gate
and drain. This input pulse drives the Lithium (Li) ions in or out of the
channel host material, WO${}_{\textbf{3}}$. Unlike conventional memristors,
these “write” (programming conductance) and “read” (measuring conductance)
processes in three-terminal electrochemical synapses are decoupled [Fig.
1(a,b)], meaning that pulses have to be directed to different terminals to
elicit the different behaviours and the read-out process can be carried out in
parallel with the write pulse. This difference, when compared to standard two-
terminal memristors where a “read” operation might affect the device
conductance, allows the use of low programming voltage and better state
retention [62]. State retention and low programming voltage are important when
designing standard memristive networks for spiking and artificial neural
network implementations as these allow for stable networks while consuming
significantly less power. Moreover, multiple small pulses can be used to
change the device conductance and retain its “state” for longer duration,
allowing the artificial networks to mimic their biological counterparts.

These three-terminal devices are specially suitable for neuromorphic networks
owing to their ability to model simple synaptic kernels such as exponential
decay and spike temporal plasticity (STP). For this reason, we focus on a
version of Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ memristor that uses a self-
gate design in which transitory effects are predominant over long-term
effects, as detailed in previous work [61]. The device behavior can be
summarised as shown in Fig. 1(b,c), where a positive ‘write’ pulse induces a
linear increase of the channel conductance (Pulse phase in blue, assuming the
electric charges accumulate linearly at the interface as a function of time)
from a starting resting phase (in green). The conductance change is
illustrated in Fig. 1(c). Positive Li+ charges driven by the write-pulse
injected electric field, accumulate at the channel/electrolyte interface and
induce charge-balancing electrons at the other side of the interface,
increasing the channel conductance. When the write pulse is removed, the
device undergoes a ‘relaxation’ phase (orange) where the charges are released
back to the electrolyte and the memristor slowly returns to its initial state
(green). The Relaxation phase can be fit using a double exponential decay
function (Eq. 19), as explained by Kohlrausch-Williams-Watts (KWW) relaxation
model [67, 68]. This behavior has also been found in other electro-chemical
devices [59] [shown in Fig. 1(b), bottom panel, relaxation phase].

This Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ device, like any analog device,
has inherent stochastic noise, which can lead to loss of information in the
encoded variables. Thus, when designing algorithms for this class of devices,
it is important to assess how this loss of precision might affect computation.
Since hardware implementations are expensive and time consuming, it is
desirable to model the device behavior and test it through simulation before
eventually moving to chip design and mass manufacturing.

###  2.1 Descriptive model

Similar to other works [67, 68, 61, 45, 59], we first develop a descriptive
model to fit the response of the device to a single square pulse stimulation.
The overall dynamics of the memristor for single pulse stimulation can be
described using Eqs. 1 and 2.

$G_{rise}(t,t_{i})=L_{1}(t,t_{i})+L_{2}(t,t_{i})$ (1)
$G_{decay}(t,t_{i})=E_{1}(t,t_{i})+E_{2}(t,t_{i})$ (2)

where $G_{rise}(t,t_{i})$ (Eq. 1) represents the conductance change of the
device channel during a positive WRITE pulse stimulus at the time $t_{i}$, and
$G_{decay}(t,t_{i})$ (Eq. 2) models the decay after the pulse is removed. The
reason of using two different expressions for both linear ($L_{1,2}$) and
exponential ($E_{1,2}$), is to keep the two conductance changes associated
with two different temporal decays independent, assuming the presence of two
separate populations of ionic charges being populated linearly at the same
rate (Eq.3), but depleted exponentially with two different temporal
coefficients $\tau_{1},\tau_{2}$ (Eq.4). This dual decay model is different
compared to our previous work [45] in which the two ionic charges were assumed
to be equal, causing the simulation to diverge from the actual device
conductance for multiple WRITE pulses (In Fig. 3 we report the RMSE of the new
and old models for the same spike train)

$\begin{cases}L_{1}(t,t_{i})=A_{1}\frac{t-t_{i}}{w}+B_{1}(t_{i})\;\;t\in
S_{rise}\\\ L_{2}(t,t_{i})=A_{2}\frac{t-t_{i}}{w}+B_{2}(t_{i})\;\;t\in
S_{rise}\\\ L_{1}(t,t_{i}),L_{2}(t,t_{i})=0\;\;t\not\in S_{rise}\\\
S_{rise}=\\{(t-t_{i})<w|t_{i-1}>t>t_{i}\\}\end{cases}$ (3)

where $t_{i}$ defines the pulse onset time and $w$ is the pulse width or
duration of the pulse. $A_{1},A_{2},\tau_{1},\tau_{2}$ are model parameters
determined experimentally, based on device dimensions and pulse parameters
including $w$ but also temperature and pulse amplitude (in Volts) [61] and
[45]. $B_{1}(t_{i})$ and $B_{2}(t_{i})$ are linear bias constants as defined
in Eq. 5.

$\begin{cases}E_{1}(t,t_{i})=(A_{1}+B_{1}(t_{i}))e^{-(\frac{t-t_{i}-w}{\tau_{1}})}\;\;t\in
S_{decay}\\\
E_{2}(t,t_{i})=(A_{2}+B_{2}(t_{i})))e^{-(\frac{t-t_{i}-w}{\tau_{2}})}\;\;t\in
S_{decay}\\\ E_{1}(t,t_{i}),E_{2}(t,t_{i})=0\;\;t\not\in S_{decay}\\\
S_{decay}=\\{(t-t_{i})\geq w|t_{i-1}>t>t_{i}\\}\ \end{cases}$ (4)

To keep track of all previous pulses before $t_{i}$, we model the “baseline”
conductance of the two populations of ions through the functions
$B_{1}(t_{i})$ and $B_{2}(t_{i})$ that store the current value of the
conductance when a new ‘WRITE’ pulse at $t_{i}$ is sent to the device.

$\begin{cases}B_{1}(t_{i})=L_{1}(t_{i},t_{i-1})+E_{1}(t_{i},t_{i-1})\\\
B_{2}(t_{i})=L_{2}(t_{i},t_{i-1})+E_{2}(t_{i},t_{i-1})\end{cases}$ (5)

Finally, when we consider the combination of the two expressions (Eq. 1,2),
during and after an $i^{th}$ pulse stimuli, we obtain $G(t)$, modelling the
total response of the channel conductance.

$G(t)=\sum_{t_{i}}G_{rise}(t,t_{i})+G_{decay}(t,t_{i})$ (6)

We fit the model to multiple recordings of pulses with different pulse
amplitude (ranging from 1$V$ to 4$V$) and duration (ranging from 200$\mu$s to
1 ms), to calculate the distribution of each parameter of the model (for
details, refer to the supplementary materials). Fig. 2 shows an example of the
descriptive model fitted to data recorded with a stimuli pulse of duration
200$~{}\mu$s and amplitude of 1$V$. The fitting results are reported in Table
1 and Table 2. The tables show the mean and standard deviations of the fitted
parameters estimates over 20 recordings for the 200us 1V pulse and 5
recordings for the remaining conditions. Values are normalized so that
baseline conductance before stimulation was zero and mean conductance after
stimulation (the sum of $A_{1},A_{2}$) was one. The decay constants,
$\tau_{1}$ and $\tau_{2}$, were found to be in the range of $10-20~{}ms$ and
around $500~{}ms$ respectively. The fitting accuracy was measured using the
root mean square error (RMSE) and referred as $\eta$. This variable, $\eta$,
is used to define the inherent noise in the device and recording setup. We
assume this noise to be additive white Gaussian noise with standard deviation
equal to the RMSE ($\eta$).

###  2.2 Generative model

To extend the descriptive model to a generative one, we consider the model
parameters $A_{1},A_{2},\tau_{1},\tau_{2}$ as stochastic variables
$\hat{A_{1}},\hat{A_{2}},\hat{\tau_{1}},\hat{\tau_{2}}$ drawn from a Gaussian
distribution with mean and standard deviation derived from the descriptive
model as shown in Eq. 7. Additionally, we add the white noise $\hat{\eta}$, to
the conductance $G_{(}t)$ at every time step $t$. For parameter $\eta$, the
variance, $\sigma_{\eta}$ is calculated as the mean deviation of the samples
from each individual fit (RMSE). Equations 3, 4, 6 can now be rewritten to
create the generative model as per Eqs. 8, 9 and 10.

$\begin{cases}\hat{A_{1}}=\mathcal{N}(\overline{A_{1}},\sigma_{A_{1}})\\\
\hat{A_{2}}=\mathcal{N}(\overline{A_{2}},\sigma_{A_{1}})\\\
\hat{\tau_{1}}=\mathcal{N}(\overline{\tau_{1}},\sigma_{\tau_{1}})\\\
\hat{\tau_{2}}=\mathcal{N}(\overline{\tau_{2}},\sigma_{\tau_{2}})\\\
\hat{\eta}=\mathcal{N}(0,\sigma_{\eta})\end{cases}$ (7)

$\begin{cases}L_{1}(t,t_{i})=\hat{A_{1}}\frac{t-t_{i}}{w}+B_{1}(t_{i})\;\;t\in
S_{rise}\\\ L_{2}(t,t_{i})=\hat{A_{2}}\frac{t-t_{i}}{w}+B_{2}(t_{i})\;\;t\in
S_{rise}\\\ L_{1}(t,t_{i}),L_{2}(t,t_{i})=0\;\;t\not\in S_{rise}\\\
S_{rise}=\\{(t-t_{i})<w|t_{i-1}>t>t_{i}\\}\end{cases}$ (8)

$\begin{cases}E_{1}(t,t_{i})=(\hat{A_{1}}+B_{1}(t_{i}))e^{-(\frac{t-t_{i}-w}{\hat{\tau_{1}}})}\;\;t\in
S_{decay}\\\
E_{2}(t,t_{i})=(\hat{A_{2}}+B_{2}(t_{i})))e^{-(\frac{t-t_{i}-w}{\hat{\tau_{2}}})}\;\;t\in
S_{decay}\\\ E_{1}(t,t_{i}),E_{2}(t,t_{i})=0\;\;t\not\in S_{decay}\\\
S_{decay}=\\{(t-t_{i})\geq w|t_{i-1}>t>t_{i}\\}\ \end{cases}$ (9)

$G(t)=\sum_{t_{i}}G_{rise}(t,t_{i})+G_{decay}(t,t_{i})+\hat{\eta}$ (10)

Figure 3 shows the comparison between the generative model (orange) and an
actual recording (light blue) for a random spike train. The difference between
the two curves is also shown (purple). The mean RMSE between the device
response and the presented model computed over eleven repetitions was
(0.413$\pm$0.027) (Unitless as the device conductance was normalized before
fitting $\pm$ standard deviation). This result shows an improvement over the
previous model presented in [45], which had a mean RMSE of 0.601$\pm$0.04 for
the same test.

##  3 Memristor as an abstraction of synaptic integration properties

Our previous work [45] proposed that these memristors can be used to model a
synapse in artificial networks due to their ability to set and change
conductance reliably, and also to respond dynamically to spiking inputs,
demonstrating biological properties like Short Time Plasticity (STP). In time
sensitive tasks (like spike patterns classification) the presented device
could boost separation of classes in single neurons up to 226x when compared
against standard (no-STP) memristors. Adding the device stochasticity into the
model reduced the gain to 114x. The proposed classification problem hinted
towards use of this device to enable (or improve) temporal integration for
spiking neural networks, but showed how the device stochasticity could impair
classification accuracy in simple networks. To make a definitive case for the
importance of these dynamical properties for synaptic devices, it is necessary
to demonstrate their use on a versatile neuromorphic architecture, and test
them on a real world problem. Moreover, this would also give us a better
insight on how inherent device stochasticity could affect more complex and
redundant architectures.

We define an “Ideal” deterministic memristor by setting the stochastic
parameters of our generative model (Eq. 10) as the mean value obtained from
the descriptive model. The ideal memristor parameters can be described in
Eq.11 as deterministic set from Eq. 7.

$\begin{cases}\hat{A_{1}}=\overline{A_{1}}\\\ \hat{A_{2}}=\overline{A_{2}}\\\
\hat{\tau_{1}}=\overline{\tau_{1}}\\\ \hat{\tau_{2}}=\overline{\tau_{2}}\\\
\hat{\eta}=0\end{cases}$ (11)

As the responses of this model are fully deterministic, it can be considered
as the baseline to compare the stochastic memristor model, allowing to isolate
the effect of noise on computation. For brevity, we refer to any network using
the deterministic model as “Ideal” network and ones using the full generative
model as “Noisy” networks.

###  3.1 HOTS

The Hierarchy of Event-Based Time-Surfaces deep temporal machine learning
alogorithm (HOTS, [52])is built on feature sets called Time-Surfaces. These
are computed using the relative timings of incoming events that are mapped to
a normalized feature space using exponential decay functions. HOTS allows
unsupervised feature extraction designed specifically for event based data,
which makes it extremely versatile for dealing with numerous neuromorphic
sensors and modalities [69, 70] and the Time-Surfaces, as previously
mentioned, are based on the interpolation of spikes with exponential decays,
similar to the output of the Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ device
(Fig. 3). Interestingly, a critical issue with the hardware implementation of
HOTS is the necessity to interpolate each event with an exponential decay
function, which is does not have an efficient digital hardware implementation.
However, the memristive dynamics of the Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$
electrochemical synapse inherently implement these decays. In this section, we
detail an implementation of HOTS trained on the N-MNIST dataset [71], which is
used as a standard in the neuromorphic community for object recognition. As
N-MNIST is a visual event based dataset, each event can be described by the
tuple:

$ev_{i}=(x_{i},y_{i},p_{i},t_{i})$ (12)

where $x_{i}$ and $y_{i}$ are the pixel positions of an event occurring.
$p_{i}$ is the polarity of the event (a bit encoding the sign of the
brightness change recorded by the pixel), and $t_{i}$ is the timestamp of the
event, usually in microseconds. This process is shown in Fig. 4 in which a
single input digit (a) is captured by an event-based sensor (b), which in turn
emits different spike train of events per pixel position (P(0,0),P(1,0),P(0,1)
etc. in (c)) when a brightness change is detected [52]. In HOTS, every
incoming event gives rise to an instantaneous rise exponential decay kernel
with no memory, whereas in our model we use a double exponential decay with
STP (we investigate further into this distinction in the supplementary
materials). Moreover, each spike is treated as a WRITE signal to a separate
synaptic device (Fig. 4d), which generates a continuous response of
conductance change (Fig. 4e). Every time a new event is captured, a spatial
neighborhood of pixels is sampled around the originating event (light blue
dashed line in Fig. 4e) and Time-Surfaces are created. This spatial
neighborhood is a squared spatial window of size $l$ encapsulating the last
event of the same polarity $p_{j}$ produced by every pixel before the
reference event $ev_{i}$. We define it by the set $\textit{T}_{i}(l,p)$:

$\begin{gathered}\textit{T}_{i}(l,p)=\max_{j\leq
i}\\{t_{j}|x_{j}=[x_{i}-l,x_{i}+l],\\\
y_{j}=[y_{i}-l,y_{i}+l],p_{j}=p\\}\end{gathered}$ (13)

As we want each Time-Surface to be fully representative of the original data,
we create a super-set $\textbf{{T}}_{i}(l)$ by concatenating all the different
$\textit{T}_{i}(l,p)$ of different polarities along an arbitrary dimension (x
or y). Finally, we can build Time-Surfaces by sampling our model $G(t)$ (Eq.
10) at the time $t_{i}$ after applying it on the set $\textbf{{T}}_{i}(l)$:

$\textit{S}_{i}(l)=G(t_{i},\textbf{{T}}_{i}(l))$ (14)

We can build this bidimensional map for any received event, and use
unsupervised learning methods as K-means, to learn local spatiotemporal
features of the input data. Now, if a new event is used to generate a Time-
Surface, we can compare it with the clusters extracted previously from the
k-means algorithm and assign the closer “centroid” to the input surface. The
index of the cluster can be used as a new polarity $p_{i}$ for the reference
event $ev_{i}$ that generated the Time-Surface, creating new sets of events
that in turn can be used to create new Time-Surfaces and new sets of features.
We can now define a layer k as a single iteration of this process so that its
input events are:

$ev^{[k]}_{i}=(x^{[k]}_{i},y^{[k]}_{i},p^{[k]}_{i},t^{[k]}_{i})$ (15)

and its output events are:

$ev^{[k+1]}_{i}=(x^{[k+1]}_{i},y^{[k+1]}_{i},p^{[k+1]}_{i},t^{[k+1]}_{i})$
(16)

Where the superscript ${[k]}$ is used to indicate the relative k layer. An
example of a 2 layer architecture is shown in Fig. 4 (from (g) to (m)). In (g)
we show the input events produced from an event-based camera and a reference
event $ev^{[1]}_{i}$. We sample the events in the neighborhood defined by
$\textbf{{T}}^{[1]}_{i}(l)$ (red square) and produce a Time-Surface
$\textit{S}^{[1]}_{i}(l)$ (h). This time surface get assigned to a cluster of
Layer 1 (i) producing a new output event $ev^{[2]}_{i}$ with a new polarity
$p^{[2]}_{i}$ (j). To increase the spatial integration, we can also subsample
the space of the new events so that:

$\begin{gathered}ev^{[2]}_{i}=(x^{[2]}_{i},y^{[2]}_{i},p^{[2]}_{i},t^{[2]}_{i})=\\\
=(x^{[1]}_{i}//u,y^{[1]}_{i}//u,p^{[2]}_{i},t^{[1]}_{i})\end{gathered}$ (17)

Where $u$ is the lateral dimension of the bigger patch aggregating the
response of multiple pixels together, and $//$ represents the “floor”
function, which outputs only the integer part of a division. The entire
process can be repeated again as shown in Fig. 4(k,l,m), up until the
preferred temporal and spatial integration is achieved.
The network cannot perform classification on its own, thus we need to pair it
with a classifier. Similarly to [52], we build a normalized histogram
$\mathcal{H}_{p}$ of the cluster responses (polarities) of the last layer to
each recording in the dataset, counting how many time events for every
polarity are produced per recording. Originally, classification was performed
by using a distance metric (Euclidean distance on absolute histogram,
Euclidean distance on normalized histograms, and Bhattacharyya) to match the
histograms between test set and train set, to find corresponding labels.
Instead, we preferred to feed the histograms as input vectors to a polynomial
Support Vector Classifier (SVC) of order 3 in order to improve the recognition
rate over the N-MNIST dataset. The classification performance is presented in
Table 3, as well as Supplementary Table 4 and 5 and further discussed in the
next section.

#### Architecture:

The architecture we propose to test is a two layered network as depicted in
Fig. 4(g to m). The Time-Surface lateral dimensions are $l^{[1]}=7$ and
$l^{[2]}=3$ respectively for the first and second layer, while the sub-
sampling parameter is $u=7$ (Eq. 17), which reduces the N-MNIST resolution
from 28x28 to 4x4. After comparing multiple simulations with different pulse
characteristics from Table 2 and 1, we chose the data generated with 1V 200us
as it performed the better over 5 consecutive runs (results shown in
Supplementary Tab. 4). During this test the cluster number was $N^{[1]}=16$
for the first layer and $N^{[2]}=512$ for the second. We then decided to run a
parameter search for the optimal number of clusters and set them to
$N^{[1]}=32$ and $N^{[2]}=64$, as it reached a good compromise between
performance and computation time (Supplementary Fig. 6 for details). The tests
were performed on 10% of N-MNIST train-set and test-set [71]. Files were
randomly selected at each run.

#### Effects of noise on network behavior:

To analyze the effects of noise on recognition rate of a full architecture in
a real word task (N-MNIST), we compare the device response in the presence of
memristor induced stochastic noise (’Noisy’ network) against an ideal non-
noisy exponential decay based network (or ‘Ideal’ network). As the
implementation of in-situ learning on the memristive chip was beyond the scope
of this paper, we limited noise analysis to inference only. Temporal feature
extraction was performed using the Ideal network model and the features were
then clustered using k-means algorithm. The Noisy network and the Ideal
network therefore share the same sets of clusters, but in the case of the
Noisy network, the Time-Surfaces generated by the test set were perturbed by
the device stochasticity we modelled using the distribution of parameters we
obtained experimentally (Eq. 7).

Thus, the performance of the Ideal and Noisy networks was analyzed using
N-MNIST test set. We found that the amount of noise in the
Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ memristor did not negatively affect
recognition rate over the average of seven randomized test runs [Fig. 5(a)].
To make sure that the Support Vector Classifier used for discrimination can
compensate for noise, we also performed a separate classification using
histogram based Euclidean distance discrimination as described by Lagorce et
al. [52] in the HOTS introduction. Thus, we created four different result
groups - Ideal and Noisy network with Support Vector Classifier (Ideal SVC and
Noisy SVC) and Ideal and Noisy network with Histogram Classifier (Ideal Hist
and Noisy Hist). Each network was trained on the same training set of N-MNIST
data and tested on randomly chosen N-MNIST test data. The classification was
also performed using either the whole network of 2 layers and a simpler single
layer network made of only the first layer. Since the single layer network
only consists of the input event layer and the memristive synapse layer, it is
directly affected by the added noise. The results (Tab. 3) show that the
addition of noise to the memristor dynamics does not significantly affect the
classification accuracy.

Even though the results did not suggest any clear effect of device
stochasticity on recognition rate, further analysis shows that the noise
affected both Time-Surface based feature generation [Fig. 5(b)] and cluster
assignment [Fig. 5(c,d,e)]. In the first case (b) we can see a Time Surface of
the whole input ($l=28$), built with using the “Ideal” model and the “Noisy”
model. The digit is still recognizable, but variation in individual pixel is
present due to “background noise” introduced because each individual device’s
conductance is exactly 0 at any point in time.

The effect of this alteration is the “dislocation”, or incorrect assignment of
Time-Surfaces (and event polarities) to the wrong cluster [Fig. 5(c,d,f)]. The
mean dislocation for a single run was around 11% of events for the first layer
and 8% of events for the second, showing that multiple layers might actually
decrease the noise effect by increasing temporal and spatial integration. Even
after this consideration, it seems hard to believe that recognition rates
would not be negatively affected by noise in the first layer [Fig. 5(a)].
Since both SVC and histogram based classifiers use the histogram of the entire
response from single recordings, they might be able to integrate enough
information to compensate for the dislocation error.

To make sure this is not the case, we calculated the overall information
provided by the spike patterns under “Ideal” and “Noisy” conditions,
irrespective of the classifier, using a method originally presented in Akolkar
et al.[72] to study how precise timings can affect spike-based pattern
recognition (More info can be found in the supplementary material). We use
this method to calculate the average Mutual Information ($MI$) between each
cluster and the input digit classes at multiple timescales and we show the
results for Ideal and Noisy architectures [Fig. 5(f-g)]. Interestingly, the
data show that architecture with the induced noise does not lose significant
information [Fig. 5(h)] at any temporal scale, suggesting that the temporal
integration performed by the histogram was not compensating the noise.
Moreover, the result suggests that event dislocation does not decrease the
mutual information, meaning that the majority of events get reassigned to
equally informative clusters.

In fact, when we look at the Mutual Information Percentage Loss [Fig. 5(h)]
(calculated as $(MI_{Ideal}-MI_{Noisy})/MI_{Ideal}$ per each individual run),
we find that for most time scales, Mutual Information loss at first layer does
not exceed 5% (for 11% of dislocated events), and 2% at second layer (for 8%
of dislocated events). In addition, when dislocated events get reassigned to
different clusters, they could increase Mutual Information of the network, as
demonstrated by the negative values of Fig. 5(h). It is important to note,
that this is true because the initial feature extraction of the Ideal network
is not perfect, allowing for some events, that were originally assigned to the
wrong cluster, to be re-assigned to the correct one, increasing the mutual
information for both original and newly assigned cluster. As we decrease the
temporal window size to a few milliseconds, the extracted Mutual Information
approaches zero, which causes dislocations to have more impact on the Mutual
Information Percentage Loss.

##  4 Discussion

Recent developments in semiconductor technology have led to the design and
creation of a new class of devices called memristors. It has been predicted
that memristors will be used in the near future as the atomic component of
more advanced and complex systems, which can provide performance superior to
conventional transistor-based hardware [73]. Deep learning architectures have
extremely high power requirements, which are often addressed by offloading
computation to distant servers. However, fully autonomous robots and drones
will require low power AI processors to operate independently. The potential
for memristors to implement neuromorphic architectures for AI has motivated a
resurgence of interest in analog computing architectures. This resurgence of
interest in neuromorphic computing has coincided with the development of new
neuromorphic sensors, such as neuromorphic skin [74] and retinas [75], which
have been evolving to fit robotic applications better.

Memristors have been used to implement components of neuromorphic
architectures. Because many memristors (RRAM [76], PCM [77], STT-MRAM [78],
FeFET [79]) do not exhibit natural temporal dynamics, they have been used
primarily for spatial integration of signals. This holds promise for
improvements in hardware implementations of neuromorphic networks, which have
been limited by the difficulty in efficiently modeling synaptic weights in
silicon. However, because they have only been used to store the weights
statically, they still require intricate circuitry to generate spike-based
temporal dynamics.

In this paper, we present a model of memristor that allows time-encoding. We
focus on a subgroup of memristors, like the proposed
Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ electrochemical memristor, which offers
biological synapse-like dynamics like exponential decay and short term
plasticity. This means that for the first time, we use the memristive device
in the computation layer to provide the dynamic temporal integration necessary
to emulate realistic synaptic processes. These will enable the implementation
of realistic synapse dynamics, and thus the advent of a new generation of
neuromorphic electronic robotic brains.

We developed a model capturing the stochasticity of the electrochemical
synapse. We used this model to implement a learning architecture based on the
work of Lagorce et al. [52], time surfaces, which are a widely used class of
event-based feature descriptors. We further demonstrated their use in HOTS, a
hierarchical network for event-based processing. The memristive dynamics were
successfully used to learn a HOTS network and perform a recognition task.

We demonstrated that the proposed method is robust to the intrinsic
stochasticity of the electrochemical device, by studying the effect of the
memristor noise on the recognition rate, the data representation, and the
amount of information. When we compared an ideal noiseless network to the
response of the modelled memristor, we found virtually no difference in the
recognition rate and mutual information (Fig. 5(a,f,g,h)), and only a minimal
effect on cluster assignment [Fig. 5(c,d,e)]. Moreover, we showed that double
exponential decays offers a computational edge over single exponential decays
as used by the original HOTS paper.

For concreteness, we have focused on one particular memristor, the
Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ electrochemical memristor. However,
similar temporal dynamic properties have been reported in other types of
memristors [23, 24, 25, 26]. This allows this work to be translated to other
existing memristive technologies.

While we have not focused on the implementation of spatial integration using
memristors here, we note that two complementary classes of memristive devices
are currently present: non-volatile and volatile. They differ in their
adequacy in spatial or temporal integration. Thus, both required elements to
move forward in the construction of electronic brains are thus available. The
advance of neuromorphic robotic brains is now tightly linked to the ability of
semiconductor manufacturing technologies to provide efficient processes able
to design hybrid systems combining transistors and memristors.

The future seems even more promising as it has been shown that memristors can
be used as a replacement of transistors. The concepts and ideas introduced in
this paper to combine both spatial and temporal integration using different
categories of memristors will make possible a new generation of neuromorphic
electronic brains, which are more compact, low power, low latency, low cost,
and widely available.

## Acknowledgments

We would like to thank Dr. Sio-Hoi Ieng for his valuable comments and
discussion during the manuscript preparation.

## References

  * [1] Harry A. Pierson and Michael S. Gashler  “Deep learning in robotics: a review of recent research”  In _Advanced Robotics_ 31.16 TaylorFrancis, 2017, pp. 821–835  DOI: 10.1080/01691864.2017.1365009
  * [2] Catherine D Schuman et al.  “A survey of neuromorphic computing and neural networks in hardware”  In _arXiv preprint arXiv:1705.06963_ , 2017
  * [3] Patrick Lichtsteiner, Christoph Posch and Tobi Delbruck  “A 128 x 128 120db 30mw asynchronous vision sensor that responds to relative intensity change”  In _2006 IEEE International Solid State Circuits Conference-Digest of Technical Papers_ , 2006, pp. 2060–2069  IEEE
  * [4] Christoph Posch, Teresa Serrano-Gotarredona, Bernabe Linares-Barranco and Tobi Delbruck  “Retinomorphic Event-Based Vision Sensors: Bioinspired Cameras With Spiking Output”  In _Proceedings of the IEEE_ 102.10, 2014, pp. 1470–1484  DOI: 10.1109/JPROC.2014.2346153
  * [5] Teresa Serrano-Gotarredona and Bernabé Linares-Barranco  “A 128$\,\times$ 128 1.5% Contrast Sensitivity 0.9% FPN 3 µs Latency 4 mW Asynchronous Frame-Free Dynamic Vision Sensor Using Transimpedance Preamplifiers”  In _IEEE Journal of Solid-State Circuits_ 48.3, 2013, pp. 827–838  DOI: 10.1109/JSSC.2012.2230553
  * [6] Guang Chen et al.  “Event-Based Neuromorphic Vision for Autonomous Driving: A Paradigm Shift for Bio-Inspired Visual Sensing and Perception”  In _IEEE Signal Processing Magazine_ 37.4, 2020, pp. 34–49  DOI: 10.1109/MSP.2020.2985815
  * [7] Giacomo Indiveri, Federico Corradi and Ning Qiao  “Neuromorphic architectures for spiking deep neural networks”  In _2015 IEEE International Electron Devices Meeting (IEDM)_ , 2015, pp. 4.2.1–4.2.4  DOI: 10.1109/IEDM.2015.7409623
  * [8] Ning Qiao et al.  “A reconfigurable on-line learning spiking neuromorphic processor comprising 256 neurons and 128K synapses”  In _Frontiers in Neuroscience_ 9, 2015, pp. 141  DOI: 10.3389/fnins.2015.00141
  * [9] Elisabetta Chicca, Fabio Stefanini, Chiara Bartolozzi and Giacomo Indiveri  “Neuromorphic Electronic Circuits for Building Autonomous Cognitive Systems”  In _Proceedings of the IEEE_ 102.9, 2014, pp. 1367–1388  DOI: 10.1109/JPROC.2014.2313954
  * [10] Mike Davies et al.  “Loihi: A Neuromorphic Manycore Processor with On-Chip Learning”  In _IEEE Micro_ 38.1, 2018, pp. 82–99  DOI: 10.1109/MM.2018.112130359
  * [11] Andreas V.. Herz, Tim Gollisch, Christian K. Machens and Dieter Jaeger  “Modeling Single-Neuron Dynamics and Computations: A Balance of Detail and Abstraction”  In _Science_ 314.5796 American Association for the Advancement of Science, 2006, pp. 80–85  DOI: 10.1126/science.1127240
  * [12] L. Abbott  “Lapicque’s introduction of the integrate-and-fire model neuron (1907)”  In _Brain Research Bulletin_ 50, 1999, pp. 303–304
  * [13] Carver Mead and Ismail Mohammed  “Analog VLSI Implementation of Neural Systems”, 1989
  * [14] Yann Lecun  “Generalization and network design strategies”  In _Connectionism in perspective_ Elsevier, 1989
  * [15] Sepp Hochreiter and Jürgen Schmidhuber  “Long Short-Term Memory”  In _Neural Computation_ 9.8, 1997, pp. 1735–1780  DOI: 10.1162/neco.1997.9.8.1735
  * [16] Andre Xian Ming Chang and Eugenio Culurciello  “Hardware accelerators for recurrent neural networks on FPGA”  In _2017 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2017, pp. 1–4  DOI: 10.1109/ISCAS.2017.8050816
  * [17] D.. Strukov et al.  “Hybrid CMOS/memristor circuits”  In _2010 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2010, pp. 1967–1970  DOI: 10.1109/ISCAS.2010.5537020
  * [18] Johannes Schemmel, Johannes Fieres and Karlheinz Meier  “Wafer-scale integration of analog neural networks”  In _2008 IEEE International Joint Conference on Neural Networks (IEEE World Congress on Computational Intelligence)_ , 2008, pp. 431–438  DOI: 10.1109/IJCNN.2008.4633828
  * [19] R. Chiu et al.  “Design and implementation of a jerk circuit using a hybrid analog–digital system”  In _Chaos, Solitons and Fractals_ 119, 2019, pp. 255–262  DOI: https://doi.org/10.1016/j.chaos.2018.12.029
  * [20] Giacomo Indiveri  “Modeling Selective Attention Using a Neuromorphic Analog VLSI Device”  In _Neural Computation_ 12.12, 2000, pp. 2857–2880  DOI: 10.1162/089976600300014755
  * [21] Jongkil Park et al.  “Hierarchical Address Event Routing for Reconfigurable Large-Scale Neuromorphic Systems”  In _IEEE Transactions on Neural Networks and Learning Systems_ 28.10, 2017, pp. 2408–2422  DOI: 10.1109/TNNLS.2016.2572164
  * [22] Daniel B. Fasnacht, Adrian M. Whatley and Giacomo Indiveri  “A serial communication infrastructure for multi-chip address event systems”  In _2008 IEEE International Symposium on Circuits and Systems_ , 2008, pp. 648–651  DOI: 10.1109/ISCAS.2008.4541501
  * [23] Takeo Ohno et al.  “Short-term plasticity and long-term potentiation mimicked in single inorganic synapses”  In _Nature materials_ 10.8 Nature Publishing Group, 2011, pp. 591–595
  * [24] Rui Yang et al.  “Synaptic plasticity and memory functions achieved in a WO3- x-based nanoionics device by using the principle of atomic switch operation”  In _Nanotechnology_ 24.38 IOP Publishing, 2013, pp. 384003
  * [25] Mohammad Taghi Sharbati et al.  “Low-Power, Electrochemically Tunable Graphene Synapses for Neuromorphic Computing”  In _Advanced Materials_ 30.36 Wiley Online Library, 2018, pp. 1802353
  * [26] Linfeng Sun et al.  “Synaptic computation enabled by joule heating of single-layered semiconductors for sound localization”  In _Nano Letters_ 18.5 ACS Publications, 2018, pp. 3229–3234
  * [27] Leon Chua, Georgios Ch Sirakoulis and Andrew Adamatzky  “Handbook of Memristor Networks”  Springer Nature, 2019
  * [28] Zhongrui Wang et al.  “Resistive switching materials for information processing”  In _Nature Reviews Materials_ 5.3 Nature Publishing Group, 2020, pp. 173–195
  * [29] Yibo Li et al.  “Review of memristor devices in neuromorphic computing: materials sciences and device challenges”  In _Journal of Physics D: Applied Physics_ 51.50 IOP Publishing, 2018, pp. 503002
  * [30] Qingzhou Wan et al.  “Emerging artificial synaptic devices for neuromorphic computing”  In _Advanced Materials Technologies_ 4.4 Wiley Online Library, 2019, pp. 1900037
  * [31] D. Kuzum, R… Jeyasingh and H.-.. Wong  “Energy efficient programming of nanoelectronic synaptic devices for large-scale implementation of associative and temporal sequence learning”  In _2011 International Electron Devices Meeting_ , 2011, pp. 30.3.1–30.3.4  DOI: 10.1109/IEDM.2011.6131643
  * [32] Manan Suri et al.  “Phase change memory as synapse for ultra-dense neuromorphic systems: Application to complex visual pattern extraction”  In _2011 International Electron Devices Meeting_ , 2011, pp. 4–4  IEEE
  * [33] S Ambrogio et al.  “Spike-timing dependent plasticity in a transistor-selected resistive switching memory”  In _Nanotechnology_ 24.38 IOP Publishing, 2013, pp. 384012  DOI: 10.1088/0957-4484/24/38/384012
  * [34] Sukru B. Eryilmaz et al.  “Brain-like associative learning using a nanoscale non-volatile phase change synaptic device array”  In _Frontiers in Neuroscience_ 8, 2014, pp. 205  DOI: 10.3389/fnins.2014.00205
  * [35] A.. Vincent et al.  “Spin-Transfer Torque Magnetic Memory as a Stochastic Memristive Synapse for Neuromorphic Systems”  In _IEEE Transactions on Biomedical Circuits and Systems_ 9.2, 2015, pp. 166–174  DOI: 10.1109/TBCAS.2015.2414423
  * [36] YeonJoo Jeong et al.  “K-means Data Clustering with Memristor Networks” PMID: 29879355  In _Nano Letters_ 18.7, 2018, pp. 4447–4453  DOI: 10.1021/acs.nanolett.8b01526
  * [37] Zhongrui Wang et al.  “Fully memristive neural networks for pattern classification with unsupervised learning”  In _Nature Electronics_ 1.2 Nature Publishing Group, 2018, pp. 137–145
  * [38] Dale Purves et al.  In _Neuroscience 2nd edition_ Sinauer Associates, 2001  URL: https://www.ncbi.nlm.nih.gov/books/NBK10855/
  * [39] Elisabetta Chicca and Giacomo Indiveri  “A recipe for creating ideal hybrid memristive-CMOS neuromorphic processing systems”  In _Applied Physics Letters_ 116.12 AIP Publishing LLC, 2020, pp. 120501
  * [40] Wulfram Gerstner, Werner M Kistler, Richard Naud and Liam Paninski  In _Neuronal dynamics: From single neurons to networks and models of cognition_ Cambridge University Press, 2014
  * [41] Michael London and Michael Häusser  “Dendritic computation”  In _Annu. Rev. Neurosci._ 28 Annual Reviews, 2005, pp. 503–532
  * [42] Yingxue Wang and Shih-Chii Liu  “Active Processing of Spatio-Temporal Input Patterns in Silicon Dendrites”  In _IEEE Transactions on Biomedical Circuits and Systems_ 7.3, 2013, pp. 307–318  DOI: 10.1109/TBCAS.2012.2199487
  * [43] Haroon Anwar, Xinping Li, Dirk Bucher and Farzan Nadim  “Functional roles of short-term synaptic plasticity with an emphasis on inhibition”  In _Current opinion in neurobiology_ 43, 2017, pp. 71–78  DOI: 10.1016/j.conb.2017.01.002
  * [44] Dean Buonomano and Tiago Carvalho  “A novel learning rule for long-term plasticity of short-term synaptic plasticity enhances temporal processing”  In _Frontiers in Integrative Neuroscience_ 5, 2011, pp. 20  DOI: 10.3389/fnint.2011.00020
  * [45] Qingzhou Wan et al.  “Low-Voltage Electrochemical LixWO3 Synapses with Temporal Dynamics for Spiking Neural Networks”  In _Advanced Intelligent Systems_ n/a.n/a, 2021, pp. 2100021  DOI: https://doi.org/10.1002/aisy.202100021
  * [46] C.. Fung, K.. Wong, He Wang and Si Wu  “Dynamical Synapses Enhance Neural Information Processing: Gracefulness, Accuracy, and Mobility” PMID: 22295986  In _Neural Computation_ 24.5, 2012, pp. 1147–1185  DOI: 10.1162/NECO“˙a“˙00269
  * [47] Luozheng Li et al.  “Dynamic Information Encoding With Dynamic Synapses in Neural Adaptation”  In _Frontiers in Computational Neuroscience_ 12, 2018, pp. 16  DOI: 10.3389/fncom.2018.00016
  * [48] Ting Chang, Sung-Hyun Jo and Wei Lu  “Short-term memory to long-term memory transition in a nanoscale memristor”  In _ACS nano_ 5.9 ACS Publications, 2011, pp. 7669–7676
  * [49] Hyungkwang Lim et al.  “Short-term memory of TiO2-based electrochemical capacitors: empirical analysis with adoption of a sliding threshold”  In _Nanotechnology_ 24.38 IOP Publishing, 2013, pp. 384005
  * [50] Radu Berdan et al.  “Emulating short-term synaptic dynamics with memristive devices”  In _Scientific reports_ 6.1 Nature Publishing Group, 2016, pp. 1–9
  * [51] Muhammad Haris, Greg Shakhnarovich and Norimichi Ukita  “Space-Time-Aware Multi-Resolution Video Enhancement”  In _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020, pp. 2856–2865  DOI: 10.1109/CVPR42600.2020.00293
  * [52] Xavier Lagorce et al.  “Hots: a hierarchy of event-based time-surfaces for pattern recognition”  In _IEEE transactions on pattern analysis and machine intelligence_ 39.7 IEEE, 2016, pp. 1346–1359
  * [53] Amos Sironi et al.  “HATS: Histograms of averaged time surfaces for robust event-based object classification”  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 1731–1740
  * [54] Antoine Grimaldi et al.  “A homeostatic gain control mechanism to improve event-driven object recognition”  In _2021 International Conference on Content-Based Multimedia Indexing (CBMI)_ , 2021, pp. 1–6  DOI: 10.1109/CBMI50038.2021.9461901
  * [55] Elliot J Fuller et al.  “Li-ion synaptic transistor for low power analog computing”  In _Advanced Materials_ 29.4 Wiley Online Library, 2017, pp. 1604310
  * [56] Yiyang Li et al.  “Low-voltage, CMOS-free synaptic memory based on Li x TiO2 redox transistors”  In _ACS applied materials & interfaces_ 11.42 ACS Publications, 2019, pp. 38982–38992
  * [57] Elliot J Fuller et al.  “Parallel programming of an ionic floating-gate memory array for scalable neuromorphic computing”  In _Science_ 364.6440 American Association for the Advancement of Science, 2019, pp. 570–574
  * [58] Chuan-Sen Yang et al.  “All-Solid-State Synaptic Transistor with Ultralow Conductance for Neuromorphic Computing”  In _Advanced Functional Materials_ 28.42 Wiley Online Library, 2018, pp. 1804170
  * [59] Jiadi Zhu et al.  “Ion gated synaptic transistors based on 2D van der Waals crystals with tunable diffusive dynamics”  In _Advanced Materials_ 30.21 Wiley Online Library, 2018, pp. 1800195
  * [60] Chuan Qian et al.  “Artificial synapses based on in-plane gate organic electrochemical transistors”  In _ACS applied materials & interfaces_ 8.39 ACS Publications, 2016, pp. 26169–26175
  * [61] Qingzhou Wan et al.  “(Bi0.2Sb0.8)2Te3 based dynamic synapses with programmable spatio-temporal dynamics”  In _APL Materials_ 7.10, 2019, pp. 101107  DOI: 10.1063/1.5106381
  * [62] Yoeri Burgt et al.  “A non-volatile organic electrochemical device as a low-voltage artificial synapse for neuromorphic computing”  In _Nature materials_ 16.4 Nature Publishing Group, 2017, pp. 414–418
  * [63] Xiahui Yao et al.  “Protonic solid-state electrochemical synapse for physical neural networks”  In _Nature communications_ 11.1 Nature Publishing Group, 2020, pp. 1–10
  * [64] Yue Li et al.  “One Transistor One Electrolyte-Gated Transistor Based Spiking Neural Network for Power-Efficient Neuromorphic Computing System”  In _Advanced Functional Materials_ Wiley Online Library, 2021, pp. 2100042
  * [65] Yiyang Li et al.  “Filament-Free Bulk Resistive Memory Enables Deterministic Analogue Switching”  In _Advanced Materials_ 32.45 Wiley Online Library, 2020, pp. 2003984
  * [66] Jianshi Tang et al.  “ECRAM as scalable synaptic cell for high-speed, low-power neuromorphic computing”  In _2018 IEEE International Electron Devices Meeting (IEDM)_ , 2018, pp. 13–1  IEEE
  * [67] Graham Williams and David C Watts  “Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function”  In _Transactions of the Faraday society_ 66 Royal Society of Chemistry, 1970, pp. 80–85
  * [68] Hua-Min Li et al.  “Electric double layer dynamics in poly (ethylene oxide) LiClO4 on graphene transistors”  In _The Journal of Physical Chemistry C_ 121.31 ACS Publications, 2017, pp. 16996–17004
  * [69] Sio-Hoi Ieng, Joao Carneiro, Marc Osswald and Ryad Benosman  “Neuromorphic Event-Based Generalized Time-Based Stereovision”  In _Frontiers in Neuroscience_ 12, 2018, pp. 442  DOI: 10.3389/fnins.2018.00442
  * [70] Gregory Gauvain et al.  “Optogenetic therapy: high spatiotemporal resolution and pattern discrimination compatible with vision restoration in non-human primates”  In _Communications biology_ 4.1 Nature Publishing Group, 2021, pp. 1–15
  * [71] Garrick Orchard, Ajinkya Jayawant, Gregory K Cohen and Nitish Thakor  “Converting static image datasets to spiking neuromorphic datasets using saccades”  In _Frontiers in neuroscience_ 9 Frontiers, 2015, pp. 437
  * [72] Himanshu Akolkar et al.  “What can neuromorphic event-driven precise timing add to spike-based pattern recognition?”  In _Neural computation_ 27.3 MIT Press One Rogers Street, Cambridge, MA 02142-1209, USA journals-info …, 2015, pp. 561–593
  * [73] Dmitri Strukov and Stan Williams  “Exponential ionic drift: Fast switching and low volatility of thin-film memristors”  In _Applied Physics A_ 94, 2009, pp. 515–519  DOI: 10.1007/s00339-008-4975-3
  * [74] Luke E. Osborn et al.  “Prosthesis with neuromorphic multilayered e-dermis perceives touch and pain”  In _Science Robotics_ 3.19 Science Robotics, 2018  DOI: 10.1126/scirobotics.aat3818
  * [75] Chiara Bartolozzi et al.  “Embedded neuromorphic vision for humanoid robots”  In _CVPR 2011 WORKSHOPS_ , 2011, pp. 129–135  DOI: 10.1109/CVPRW.2011.5981834
  * [76] Wei Wang et al.  “Learning of spatiotemporal patterns in a spiking neural network with resistive switching synapses”  In _Science advances_ 4.9 American Association for the Advancement of Science, 2018, pp. eaat4752
  * [77] Abu Sebastian et al.  “Tutorial: Brain-inspired computing using phase-change memory devices”  In _Journal of Applied Physics_ 124.11 AIP Publishing LLC, 2018, pp. 111101
  * [78] Julie Grollier, Damien Querlioz and Mark D Stiles  “Spintronic nanodevices for bioinspired computing”  In _Proceedings of the IEEE_ 104.10 IEEE, 2016, pp. 2024–2039
  * [79] Sourav Dutta et al.  “Supervised Learning in All FeFET-Based Spiking Neural Network: Opportunities and Challenges”  In _Frontiers in Neuroscience_ 14 Frontiers, 2020, pp. 634
  * [80] Vandana Padala, Arindam Basu and Garrick Orchard  “A noise filtering algorithm for event-based asynchronous change detection image sensors on truenorth and its implementation on truenorth”  In _Frontiers in neuroscience_ 12 Frontiers, 2018, pp. 118
  * [81] Adam Coates and Andrew Y. Ng  “Learning Feature Representations with K-Means”  In _Neural Networks: Tricks of the Trade: Second Edition_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 561–580  DOI: 10.1007/978-3-642-35289-8˙30
  * [82] David Sculley  “Web-scale k-means clustering”  In _Proceedings of the 19th international conference on World wide web_ , 2010, pp. 1177–1178

## Supplementary materials

### Model fitting and parameters estimation

After defining the model, we need to estimate its parameters, namely,
$A_{1},A_{2},\tau_{1},\tau_{2}$. We apply a single voltage pulse at the Gate
using a function generator (Tektronix AFG3252C). We then read the channel
current IDS by applying a small reading voltage bias (100 mV) between source
(S) and drain (D) using a Keithley Semiconductor Parameter Analyzer (4200-SCS)
at a sampling rate of 6 ms. We calculate the device conductance using Ohm’s
rule; GDS = IDS/VDS. This temporal precision is not enough to evaluate
$G_{rise}$, as all the pulses duration are sub-millisecond. Therefore we can
only determine the parameters by fitting $G_{decay}$. Since we are using a
single WRITE pulse per each recording the previous state of the device can be
simplified $B_{1}(t_{0})+B_{2}(t_{0})=0$. The device conductance is now
defined by:

$G_{rise}(t,t_{0})=(A_{1}+A_{2})\frac{t-t_{0}}{w}$ (18)
$G_{decay}(t,t_{0})=A_{1}e^{-(\frac{t-t_{0}-w}{\tau_{1}})}+A_{2}e^{-(\frac{t-t_{0}-w}{\tau_{2}})}$
(19)

$G(t)=\begin{cases}G_{rise}(t,t_{0})\;\;when\;\;(t-t_{0})<w|t>t_{0}\\\
G_{decay}(t,t_{0})\;\;when\;\;(t-t_{0})\geq w|t>t_{0}\\\
0\;\;otherwise\end{cases}$ (20)

Before fitting, we scale the data by removing the mean base conductance at the
steady state, and we normalize by dividing the data by the mean maximum
conductance over multiple recordings. We took this choice to simplify the data
representation both during analysis and simulation of the network that will be
addressed later in the paper. After this step, we fit $G_{decay}$ using the
Gauss-Newton algorithm in the statistical analysis software JMP®. We can now
obtain sets of all the necessary parameters from multiple recordings of the
$Li_{x}WO_{3}$ synaptic memristor.

### Double exponential dynamics vs single exponential dynamics

Models like HOTS were designed with a no memory, single decaying exponential
dynamic, making each layer sensitive only to only one temporal scale at a
time. For this reason, to integrate temporal information, HOTS architectures
need multiple layers with increasingly longer time constants [52]. The
proposed electrochemical device, however, elicits a double exponential decay
response, that could potentially extract temporal information at multiple
timescales at once. When we build a two layered HOTS model, with the same
decays as generated by the electrochemical device for an input pulse of
200$\mu$s and 1V (Layer 1 $\tau_{1}=5ms$ and Layer 2 $\tau_{1}=92ms$) we find
comparable results on the N-MNIST dataset (Tab. 5).

Overall, the Ideal memristor performs better, and we believe for two reasons:
1) the ability of double exponential dynamics to extract temporal information
at multiple scales; 2) the effect of Short Term Plasticity in reducing noise.
To study these two effects, we compared Time-Surfaces between HOTS and the
Ideal Memristor Model Fig. 8. In $(a)$ we show the full 28x28 Time-Surface of
a single digit represented with : a single exponential decay with no memory
(for $\tau=5ms$ and $\tau=92ms$) and a double exponential decay + STP (Ideal
memristor model). In $(b)$ we show the same 11x11 portion for all three
examples, in which we identify three different regions in shades of blue: the
region of recent events (light blue), the region of sensor noise (blue) and
the region of previous events (dark blue). The noise described in this section
is not related to the memristor dynamics instead it is a form of salt and
pepper noise, common in event based cameras [80]. We use standard deviation
(SD) $(c)$ as a metric of how much of information (pixel value variation) or
noise is present in the selected regions. When compared to single exponential
decay, the double exponential decay in addition to the STP reduces the gap
between recent and previous event regions therefore integrating information at
multiple timescales, leading to reduction in the effect of sensor-induced
noise.

### Mutual Information

Information was calculated using this equation:

$I_{\delta}(S;R)=\sum_{s}P(s)\sum_{r}P(r|s)_{\delta}log_{2}(\frac{P(r|s)_{\delta}}{P(r)_{\delta}})$
(21)

with $r$ being the response of the network, either a 1 or a 0, representing
the presence of an event $ev^{[k]}_{i}$ with a given polarity $p$, and $s$
being the “type” of the stimuli eliciting the response (the N-MNIST digit). In
other words, $I_{\delta}(S;R)$ gives a measure of how reliably a single
cluster of the network, tends to respond to the same type of input stimuli.
Since we used an equal number of different digits to calculate $I_{\delta}$
the probability of $s$ is uniform $P(s)=1/10$. For $P(r)_{\delta}$, we
randomly extracted events $ev^{[k]}_{i}$ from single recordings, and then
calculated the probability of finding an event $ev^{[k]}_{j}$ with the same
polarity, in any other recording, within a time window of length $\delta$
centered around the time stamp $t_{i}$ of the original event. Similarly to
Akolkar et al. [72] responses were binary, so if multiple events with the same
polarity were present within the time window of the same recording, only one
was considered. To calculate $P(r|s)_{\delta}$ we repeated the process used
for $P(r)_{\delta}$ but only searched events over multiple recordings of the
same digit $s$. Finally, we calculated $I_{\delta}(S;R)$ over different time
windows sizes ($\delta$) and averaged the response for all the different
polarities and multiple runs, defining a metric ($MI$) of how much information
can be extracted from the network response in a given time window [Fig.
5(f,g)]. The results we show were calculated on the same recordings used for
the individual runs in Fig. 5(a), in order to offer a true comparison between
the recognition rate and mutual information. In order to validate the use of
the proposed Mutual Information calculation as a proper metric for the
stochasticity of the device and its effects on recognition rates, we simulated
an hypothetical device with the same parameters used for the results shown in
Figure 5, but increasing the standard deviation of the model parameters by 2x,
5x, 10x. These results are reported in Figure 7 and show that both Mutual
Information and recognition rate decrease inversely to higher degree of
variation in the device response (higher standard deviation in the stochastic
parameters Eq. 7).

### Short Term Plasticity and Time-Surface Normalization

While we make the argument that STP could play a role in averaging and
removing noise, it also presents a problem for clustering algorithms as
k-means. in HOTS like algorithms, Time-Surfaces are naturally bounded between
0 and 1; however, in this model, the memristor conductance is the summation of
multiple events, therefore a high rate of spiking events for the same pixel
could drive the Time-Surface representation to high values, this could impair
classification, driving the cluster algorithm to identify different clusters
for magnitude variation of the same object [81]. To avoid this problem, we
decided to normalize Time-Surfaces before clustering with k-means.

### Additional methods

All the simulation were performed in Python. Clustering was achieved using a
mini-batch k-means algorithm as a faster and efficient Python multi-threaded
implementation [82].

Figure 1: The schematic of the physical Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$
electrochemical synapse used for this paper (a) and a graphical representation
of the short-term plasticity response [(b),(c)]. (a) The
Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$ electrochemical synapse has a three-
terminal [(S)ource, (G)round, (D)rain] structure. A WRITE voltage pulse is
applied between the (G) and (D), while a small DC READ voltage (0.1 V) is
applied between the (S) and (D) to monitor the channel conductance. The
Reservoir (both Li ions and electrons) and the Channel [in light red (a),(c)]
are made of Li${}_{\textbf{x}}$WO${}_{\textbf{3}}$, which enables the exchange
flow of Li+ ions between gate and channel through the electrolyte gel.
Whenever the device receives with a WRITE voltage pulse, a transitory
conductance change is elicited, as shown in the lower graph of Fig. (b). The
resulting short term effect is dominated by an ionic gating (EDL) effect as
depicted in the lower three panels (c). Whenever a WRITE voltage Pulse is
injected, the resting equilibrium of charges [Rest phase in green (b),(c)] is
disrupted, causing a fast and linear increase in conductance [Pulse phase in
blue (b),(c)]. When the WRITE voltage stimulation is removed, the channel
conductance relaxes back with a characteristic double exponential decay
($A_{1}e^{{-t}/{\tau_{1}}}+A_{2}e^{{-t}/{\tau_{2}}})$ [Relaxation phase in
yellow (b),(c)]. Figure 2: Fitting of a single pulse experimental data
recorded from a sample device with the proposed linear rise two exponential
decay model. We scale data points so that the mean peak conductance of
multiple recordings of the same pulse is equal to 1. The rest conductance (in
green) is defined as the mean value of the points before pulse onset $t_{0}$.
To simplify fitting, it is considered to be zero. The model parameters
($A_{1},A_{2},\tau_{1},\tau_{2}$) are obtained by extracting the data after
the end of the squared pulse stimuli ($t_{0}+w$) and fitting $G_{decay}(t)$
(in yellow) using the Gauss-Newton algorithm. By knowing the length of the
pulse stimuli $w$ and the fit parameters, we can infer $G_{rise}$ (in blue)
and model the entire device response. The data used for this plot is obtained
with a 1V pulse and $w=$200us. The fitted model parameters are reported in
Tab. 1,2. Figure 3: Our generative model (orange) compared to a recording of
the device (blue) stimulated with the same spike sequence (light blue).
Differences are mainly due to the intrinsic stochasticity of our model and the
memristic device, and are prevalent in the instantaneous response of the
device to spike stimulation (as shown by the error in purple) . The Root Mean
Square Error (RMSE over 11 different comparisons was 0.413$\pm$0.27 (Unitless
$\pm$ standard deviation). This result shows an improvement over the previous
model presented in [45], which presents a mean RMSE of 0.601$\pm$0.04 for the
same test.  Table 1: Model parameters obtained with a fixed 1V pulse amplitude
and varying pulse width.

1V | $\overline{A_{1}}$ | $\overline{\tau_{1}}$ | $\overline{A_{2}}$ | $\overline{\tau_{2}}$ |  $\eta$ (RMSE)
---|---|---|---|---|---
200us | 0.57$\pm$.27 | 5ms$\pm$2ms | 0.5$\pm$.05 | 92ms$\pm$18ms | 0.11
500us | 0.74$\pm$.18 | 16ms$\pm$9ms | 0.26$\pm$.06 | 588ms$\pm$31ms | 0.12
750us | 0.78$\pm$.19 | 10ms$\pm$2ms | 0.23$\pm$.03 | 513ms$\pm$98ms | 0.04
1ms | 0.75$\pm$.29 | 10ms$\pm$3ms | 0.22$\pm$.02 | 390ms$\pm$68ms | 0.02

Each parameter is represented with by its mean value$\pm$ Standard Deviation.
$\eta$ is represented by the mean standard deviation of the individual fits,
which can also serve as a measure of the fitting precision as Root Mean Square
Error (RMSE)

Table 2: Model parameters obtained with a fixed 200us pulse width and varying
pulse amplitude.

200us | $\overline{A_{1}}$ | $\overline{\tau_{1}}$ | $\overline{A_{2}}$ | $\overline{\tau_{2}}$ |  $\eta$ (RMSE)
---|---|---|---|---|---
1V | 0.57$\pm$.27 | 5ms$\pm$2ms | 0.5$\pm$.05 | 92ms$\pm$18ms | 0.11
2V | 0.54$\pm$.29 | 7ms$\pm$1ms | 0.35$\pm$.02 | 122ms$\pm$20ms | 0.06
3V | 0.77$\pm$.24 | 13ms$\pm$6ms | 0.23$\pm$.02 | 373ms$\pm$98ms | 0.07
4V | 0.75$\pm$.17 | 11ms$\pm$3ms | 0.25$\pm$.01 | 501ms$\pm$101ms | 0.05

Each parameter is represented with by its mean value$\pm$ Standard Deviation.
$\eta$ is represented by the mean standard deviation of the individual fits,
which can also serve as a measure of the fitting precision. as Root Mean
Square Error (RMSE)

Figure 4: Mapping spatio-temporal event patterns into time surface features
using memresistive synapses (Fig. a - f), and N-MNIST classification using
full HOTS network with memristors (Fig. g - m). An input digit (a) is
presented to the event-based sensor (b), that produces asynchronous event
streams at each individual pixel, based on the luminance variation over time.
Each pixel stream is the input for a memristor (d), which interpolates the
spike train with an exponential decay kernel (in red) defined by the temporal
coefficients $\tau_{1},\tau_{2}$. Every time a new event $ev_{i}$ is sent to
any synaptic device, a new time surface is sampled, by taking the
instantaneous values across multiple pixels (in light blue figure e). The
result is a 2D map that encodes temporal correlations between events at
different pixels (f). The neighborhood of pixels around a reference event can
also be seen in (g), defined as the area delimited by the red square. We
produce this figure by coloring each pixel with the polarity of the last event
within a 10ms time window. The original events produced by the sensor have two
polarities, to indicate increase or decrease of luminance. This categorical
information represents the output of each layer of our network. Once we
produce a time surface (h) we find the closest cluster (i) and produce a new
event at the same pixel. The output event has the same timestamp as the input
event, and polarity determined by the index of the closest cluster (j). This
new set of events is spatially sub-sampled (Eq. 16) and used to generate new
time surfaces (k) that are again classified by the nearest cluster (l,m). This
process can be repeated multiple times to increase the temporal and spacial
scale of the classified features. The clusters are generated using the k-means
clustering algorithm on a training set of time surfaces.  Figure 5:
Recognition rates of the Ideal network against Noisy, with a Support Vector
Classifier SCV, and Euclidean Histogram classifier. Average results ($\pm$
Standard deviation) over seven runs are reported in Tab. 3. In (b) we show a
qualitative comparison between two Time-Surfaces of the whole input. The first
one built with the Ideal model while the second with the Noisy one. The device
inherent stochasticity affects both values of active pixels as well as
inactive background pixels, while the digit remains visible. In (c,d,e) we
show the effect of Time-Surface perturbation on network operation. When
comparing Ideal and Noisy responses to the same recording, most events (blue)
keep the assigned polarity, whereas others lose the Ideal polarity (orange) to
get assigned to a new one (green), we call this effect Cluster dislocation
(e). We analyze the Mutual Information between cluster response and N-MNIST
digits for layer 1 (f) and layer 2 (g) at different Temporal Integration
scales, to show that Cluster Dislocation do not affect network information on
average. In fact, when we look on the loss of information on individual runs
(h), we discover that event dislocation effect is small (except when temporal
integration and mutual information approach zero), and at this level of noise
Cluster Dislocation has similar chance to decrease Mutual Information as to
increase it, causing the Mutual information Loss to be negative (h)). Table 3:
Classification accuracy of Ideal network and Noisy Network with two different
classifiers (Hist,SVC)

| Ideal Hist | Noisy Hist | Ideal SVC | Noisy SVC
---|---|---|---|---
First Layer | $61.40\%\pm 1.75\%$ | $61.30\%\pm 1.71\%$ | $72.81\%\pm 4.07\%$ | $71.80\%\pm 4.69\%$
Whole Net | $76.40\%\pm 1.95\%$ | $76.44\%\pm 1.80\%$ | $87.61\%\pm 1.12\%$ | $87.76\%\pm 1.17\%$

Results obtained with $N^{[1]}=32$ clusters for layer 1 and $N^{[2]}=64$ for
layer 2 over 7 consecutive runs.

Table 4: N-MNIST classification accuracy with memristor dynamics using
different pulse parameters

| 200$\mu$s | 500$\mu$s | 750$\mu$s | 1ms
---|---|---|---|---
1V | $\mathbf{91.80\%\pm 0.88\%}$ | $90.72\%\pm 1.02\%$ | $91.73\%\pm 1.37\%$ | $90.89\%\pm 1.32\%$

Results with pulse characteristics from Table. 1

| 1V | 2V | 3V | 4V
---|---|---|---|---
200$\mu$s | $\mathbf{91.80\%\pm 0.88\%}$ | $91.21\%\pm 1.33\%$ | $90.92\%\pm 1.62\%$ | $91.03\%\pm 1.5\%$

Results obtained with $N^{[1]}=16$ clusters for layer 1 and $N^{[2]}=512$ for
layer 2 over 5 consecutive runs and with pulse characteristics from Table. 2

Figure 6: In this figure we show the results N-MNIST accuracy by number of
clusters, for a single layer (circles), and a second layer architecture
(triangles). Layer 2 results were obtained using 32 clusters in layer 1. As
recognition rate quickly reaches an asymptote around 90% recognition rate we
decided that an architecture with 32 clusters for layer 1 and 64 clusters for
layer 2 was a good compromise between performance and computation time.
Figure 7: Mutual Information and recognition rates for hypothetical memristors
stimulated at 1v 200us but with 2x, 5x, 10x the standard deviations of the
stochastic parameters. The plots show that higher degrees of stochasticity can
easily affect Mutual Information of the network (a-b) and Recognition rates
(c) for both the tested classifiers (SVC and Hist). Table 5: Classification
accuracy using memristors dynamics vs original exponential decay kernel of
HOTS

| HOTS | Ideal Mem
---|---|---
First Layer | $47.66\%\pm 8.75\%$ | $74.13\%\pm 2.61\%$
Whole Net | $89.98\%\pm 1.29\%$ | $91.80\%\pm 0.88\%$

Results obtained with $N^{[1]}=16$ clusters for layer 1 and $N^{[2]}=512$ for
layer 2 over 5 consecutive runs.

Figure 8: Figure showing the effect of multiple exponential decays and memory
on the Time-Surface representation. In (a) we show the same digit represented
with a Single exp decay with no memory (HOTS) at 5ms and 92ms, and a double
exponential decay with STP and the same temporal decays (Ideal Memristor). In
(b) we extract a 11x11 tile from (a) and identify three region of interest:
recent events (light blue), previous events (dark blue) and sensor noise
(blue). We use Standard Deviation as a measure of information and noise. A
single fast decay of 5ms can represents the temporal patterns of the recent
activity (light blue) but quickly forgets older events (dark blue). Conversely
a slower decay of 92ms retains temporal information on a longer scale (dark
blue) but relations on a smaller timescale (light blue) are lost. The
electrochemical device is capable of representing temporal information at
multiple timescales with the same pulse characteristics, while also limiting
noise caused by single random events thanks to the STP averaging multiple
events during time.


FACTORIZATION OF PLATONIC POLYTOPES INTO CANONICAL SPHERES 111appeared in
Geombinatorics, XXXI(1)(2021) 5–9.

Richard H. Hammack
Department of Mathematics and Applied Mathematics
Virginia Commonwealth University
Richmond, VA 23284-2014 USA
<EMAIL_ADDRESS>
Paul C. Kainen
Department of Mathematics and Statistics
Georgetown University
Washington, D.C. 20057 USA
<EMAIL_ADDRESS>
###### Abstract

Factorization into spheres is achieved for skeleta of the simplex, cube, and
cross-polytope, both explicitly and using Keevash’s proof of existence of
designs.

Key Phrases: Decomposition of skeleta, Steiner triples, Hanani quadruples

Decomposing a polytope’s $2$-skeleton (i.e., partitioning its $2$-cells) into
closed manifolds was studied in [1], [13], [5], [7], [9], while in [6] we
found factorizations of hypercube $2$-skeleta into boundaries of (pairwise
isomorphic) $3$-cubes. Here we obtain such sphere factorizations for
$k$-skeleta of all (non-exceptional) Platonic polytopes both explicitly for
low values of $k$ and, as a consequence of Keevash’s result [11], also
existentially with $k\geq 1$ arbitrary. Note we only use the case of
multiplicity $\lambda=1$ of [11].

Let $\Delta_{n}$ denote the $n$-dimensional simplex, whose 1-skeleton is the
graph $K_{n+1}$. Let ${\cal O}_{n}$ denote the $n$-dimensional cross-polytope,
whose 1-skeleton is the graph $K_{2n}-F$, where $F$ is a 1-factor. Let ${\cal
Q}_{n}$ denote the $n$-dimensional hypercube, whose 1-skeleton is the graph
$Q_{n}$, the $n$-fold graph Cartesian product of $K_{2}$. These are the three
families of Platonic polytopes.

For any $n$-dimensional polytope $A$ and nonnegative integer $k$, let $A^{k}$
denote the $k$-skeleton of $A$; see, e.g., [2]. If $L$ and $K$ are polytopal
$\ell$-complexes, let $L|K$ mean that L factors K; that is, $K$ and $L$ are of
equal dimension $\ell$ and $K$ is the union of subcomplexes
$L_{1},\ldots,L_{r}$, each isomorphic to $L$, such that every $\ell$-face of
$K$ is contained in exactly one of the $L_{i}$. An $\ell$-complex is even if
each $(\ell{-}1)$-face is in a positive even number of $\ell$-faces.

For $X\in\\{\Delta,{\cal O},{\cal Q}\\}$ we study the factorization of
$K=X_{n}^{\ell}$ by $L=X_{\ell+1}^{\ell}$ for $\ell\leq n{-}1$. We call $L$
the canonical $\ell$-sphere as it is the boundary of a type $X$ polytope of
dimension $\ell{+}1$. For $\ell\geq 1$, such factorizations always exist for
${\cal O}$ as the proof of the first theorem below shows explicitly.

Observe, however, that these canonical spheres might not be composed of the
minimum number of $\ell$-faces. Indeed, by Theorem 1 below, the 1-skeleton of
${\cal O}_{3}$ is factored into three 4-cycles, rather than four 3-cycles.

Clearly, $X_{n}^{\ell}$ even is a necessary condition for a factorization into
spheres. The $\ell$-skeleton of an $n$-simplex (or an $n$-cube) is an even
complex if and only if $n{-}\ell{+}1$ (the number of $\ell$-faces containing
an $(\ell{-}1)$-face) is positive and even; for the cross-polytope ${\cal
O}_{n}$, the $\ell$-skeleton is even for all $1\leq\ell<n$. Evenness is
automatic for ${\cal O}_{n}^{\ell}$ as each $(\ell{-}1)$-face is a simplex and
is part of two $\ell$-faces (also simplexes) one for each vertex in a copy of
$\overline{K_{2}}$.

The following extends [5] for $\ell=2$; cf. Spreer [13], which partitions
${\cal O}_{n}^{2}$ into face-disjoint surfaces that contain the 1-skeleton
${\cal O}_{n}^{1}$; such surfaces are called 1-Hamiltonian. Observe that our
result does not depend on [11].

######  Theorem 1.

If $1\leq\ell<n$, then ${\cal O}_{\ell+1}^{\ell}\,|\,{\cal O}_{n}^{\ell}$.

###### Proof.

For each $k\geq 1$, $\partial{\cal O}_{k}$ is the iterated topological join of
$k$ copies of the 2-point complex $\overline{K_{2}}$, which is a
$(k{-}1)$-sphere. The $\ell$-skeleton of ${\cal O}_{n}$ is the $\ell$-cell-
disjoint union of the ${{n}\choose{\ell+1}}$ copies of ${\cal
O}_{\ell+1}^{\ell}$ formed by the iterated join of each $(\ell{+}1)$-element
subset of the set $n\,\overline{K_{2}}$. ∎

For the simplex, we use the theory of Steiner-type configurations
$(v,k,\ell)$, which consist of a $v$-element set $X$ endowed with a family of
$k$-element subsets of $X$ (called blocks) such that each $\ell$-subset of $X$
is in exactly one block. The case $k=\ell{+}1$ will be what we apply to factor
$\Delta_{v-1}^{\ell-1}$ by $\Delta_{\ell}^{\ell-1}$.

The translation from design to combinatorial factorization works for special
cases, such as the Steiner triple systems as we noted in [6]. But by Keevash
[11] (see also Gowers [4] and Kalai [10]), this sort of factorization must
occur almost always for the $\ell$-skeleta of simplexes and, hence, of cubes,
as we shall explain. Write $a\,|\,b$ to mean $b$ is an integer multiple of
$a$.

Define the divisibility set (e.g., [8]) to be the set of all feasible values

${\cal
D}(k,\ell):=\Big{\\{}v:{{k-h}\choose{\ell-h}}\,\Big{|}\,{{v-h}\choose{\ell-h}},\,0\leq
h\leq\ell{-}1<k\leq v\Big{\\}},$ (1)

By [11], for $1\leq\ell\leq k$, $\exists$ finite set ${\cal
X}(k,\ell)\subset{\cal D}(k,\ell)$ such that

$\exists(v,k,\ell)\mbox{{\rm-configuration}}\iff v\in{\cal
D}(k,\ell)\setminus{\cal X}(k,\ell).$ (2)

$0$$1$$2$$3$$4$$5$$6$$7$$0$$1$$2$$3$$4$$5$$6$$7$$0$$1$$2$$3$$4$$5$$6$$7$$0$$1$$2$$3$$4$$5$$6$$7$$0$$1$$2$$3$$1$$2$$4$$7$$1$$2$$5$$6$$0$$2$$4$$6$
Four rotations of this configuration by angles of $k\frac{\pi}{2}$ for
$k=0,1,2,3$. $0234$ 2345
4567
6701 Four rotations of this configuration by angles of $k\frac{\pi}{2}$ for
$k=0,1,2,3$. 1247
3461
5603
7025 Four rotations of this configuration by angles of $k\frac{\pi}{4}$ for
$k=0,1,2,3$. 1256
2367
3470
4501 Two rotations of this configuration by angles of $0$ and $\frac{\pi}{4}$
0246
1357 Figure 1: Factoring $\Delta_{7}^{2}$ into 14 $2$-face-disjoint
tetrahedral boundaries $\Delta_{3}^{2}$

We interpret the combinatorics geometrically as applying to faces of the
simplex. Since the $(j{+}1)$-sets from the $n{+}1$ vertices of $\Delta_{n}$
correspond to the $j$-faces of $\Delta_{n}$, we obtain a factorization into
canonical spheres.

######  Theorem 2.

For $1\leq\ell$ and, assuming Keevash’s result [11],

$\Delta_{\ell{+}1}^{\ell}\,\Big{|}\,\Delta_{n}^{\ell}\,\iff\,n+1\,\in\,{\cal
D}({\ell{+}2},\ell{+}1)\setminus{\cal X}({\ell{+}2},\ell{+}1).$ (3)

Any $(v,3,2)$-configuration is a Steiner triple system (as observed in [11]),
and corresponds to a factorization of the $1$-skeleton of $\Delta_{v-1}$
(i.e., $K_{v}$) into edge-disjoint 3-cycles, the canonical (and minimum)
spheres in 1-dimension.

For the case $(v,4,3)$, corresponding to the $(v{-}1)$-simplex, Hanani
constructed a family of $4$-sets (quadruples) that cover each $3$-set exactly
once, and he proved in [8] that they exist for $v\geq 4$ if and only if
$v\equiv 2$ or $4$ (mod-$6$), so again independent of  [11], we can
characterize when the $2$-skeleton of an $n$-simplex is factored by the
boundary of a tetrahedron:

$\Delta_{3}^{2}\,|\,\Delta_{n}^{2}\;\iff\;n\equiv
1\,\mbox{or}\;3\;(\mbox{mod}\;6),\,n\geq 3.$ (4)

See Fig. 1 for an illustration of the first non-trivial case of (4).

In order to factor cube-skeleta, we use an “exponentiation” method that goes
back to Danzer [3], [12]. If $\sigma=(i_{1},\ldots,i_{k+1})$ is any
$k$-simplex in the $n$-simplex with vertex set
$[n{+}1]:=\\{1,\ldots,n{+}1\\}$, let

$2^{\sigma}:=A_{1}\times\cdots\times A_{n+1},$ (5)

where $A_{i}=[0,1]$ if $i\in\sigma$, $A_{i}=\\{0,1\\}$ if $i\notin\sigma$, and
$\times$ is topological product. Thus, $2^{\sigma}$ is a family of $2^{n-k}$
pairwise-vertex-disjoint $(k{+}1)$-dimensional cubes in the $(n{+}1)$-cube.
Let $2^{K}$ be the union of the $2^{\sigma}$ over all simplexes $\sigma\in K$
for $K$ a simplicial complex. Then we have (see Kühnel & Schulz [12])

$\mbox{{\rm(a)}}\;\;2^{\Delta_{n-1}^{\ell-1}}\,=\,{\cal
Q}_{n}^{\ell};\;\;\;\;\mbox{{\rm(b)}}\;\;2^{\partial
K}\,=\,\partial\,(2^{K});\;\;\mbox{{\it and}}\;\;$ (6)

$2^{K}\;\mbox{{\it is a manifold}}\iff K\;\mbox{{\it is a combinatorial
sphere.}}$ (7)

By (6), (a) and (b), a factorization of the $(\ell{-}1)$-skeleton of the
$(v{-}1)$-simplex by the boundary of $\Delta_{\ell}$ maps under the
exponential correspondence to a factorization of the $\ell$-skeleton of the
$n$-cube into subcomplexes isomorphic to $\partial({\cal Q}_{\ell+1})$. We
used this in [6] to factor the 2-skeleton of a cube by the boundary of a
3-cube exactly when Steiner triples exist. In [7] we used (7) to show that,
for $n$ even, the Hamiltonian factorization of the graph
$K_{n+1}=\Delta_{n}^{1}$ into spanning 1-spheres yields a factorization of the
2-skeleton of ${\cal Q}_{n+1}$ into orientable minimum-genus surfaces which
are 1-Hamiltonian.

Similarly, Hanani’s factorization (4) of the 2-skeleton of the simplex gives

${\cal Q}_{4}^{3}\,|\,{\cal Q}_{n}^{3}\iff n\geq 4,\;n\equiv
2\;\mbox{or}\;4\;(\mbox{mod}\;6).$ (8)

Again, (8) doesn’t depend on [11], while by Theorem 3, the exponential
correspondence (5), and using [11], we have

######  Theorem 3.

For $1\leq\ell$, if $n\in{\cal D}(\ell{+}1,\ell)\setminus{\cal
X}(\ell{+}1,\ell)$, then $\,{\cal Q}_{\ell+1}^{\ell}\,|\,{\cal Q}_{n}^{\ell}$.

We conjecture that the condition is also necessary. Note that the result for
$\ell=1$ follows from the true though not stated case of Theorem 3 for
$\ell=0$ since ${\cal D}(2,1)$ is the set of positive even integers and ${\cal
X}(2,1)=\emptyset$.

By the results of [11], factorability into canonical spheres is ubiquitous for
skeleta of Platonic polytopes. A geometric explanation for this phenomenon (in
special cases) would imply the existence theorem (in those cases).

## References

  * [1] U. Betke, Ch. Schulz, & J. M. Wills, Zur Zerlegbarkeit von Skeletten, Geometriae Dedicata, 5 (4) (1976) 435–451.
  * [2] H. S. M. Coxeter, Regular Polytopes, Dover Publ., New York, 1973.
  * [3] L. Danzer, Regular incidence-complexes and dimensionally unbounded sequences of such, I., Ann. of Discr. Math. 20 (1984) 115–127.
  * [4] W. T. Gowers, Probabilistic combinatorics and the recent work of Peter Keevash, Bull. Amer. Math. Soc. 54 (2017) 107–116.
  * [5] R. H. Hammack & P. C. Kainen, On $2$-skeleta of Platonic polytopes, Bull. Hellenic Math. Soc. 62 (2018) 94–102.
  * [6] R. H. Hammack & P. C. Kainen, Sphere-decompositions of hypercubes, Art of Discrete Appl. Math. 3 (2020) #P2.09.
  * [7] R. H. Hammack & P. C. Kainen, A new view of hypercube genus, American Math. Monthly 128(4)(2021) 352–359.
  * [8] H. Hanani, The existence and construction of balanced incomplete block designs, Ann. Math. Statist. 32 (1961) (2) 361–386.
  * [9] P. C. Kainen, On $2$-skeleta of hypercubes, Art of Discrete Appl. Math. 3 (2020) #P2.06
  * [10] G. Kalai, Designs exist! [after Peter Keevash], Seminaire Bourbaki 67eme annee, 2014–2015, no. 1100, 2015.
http://www.bourbaki.ens.fr/ TEXTES/ 1100.pdf

  * [11] P. Keevash, On the existence conjecture, arXiv 1401.3665v3, Aug. 5, 2019. https://arxiv.org/pdf/1401.3665.pdf
  * [12] W. Kühnel & Ch. Schulz, Submanifolds of the cube, Applied Geometry and Discrete Mathematics (The Victor Klee Festschrift), DIMAC Series in Discrete Math. and Theor. Computer Sci., 4 (1991) 423–432.
  * [13] J. Spreer, Partitioning the triangles of the cross polytope into surfaces, Beitr. Algebra Geom. 53 (2012) 473–486.