Event-Driven Contrastive Divergence for Spiking Neuromorphic Systems
Emre Neftci ∗1 , Srinjoy Das1,2 , Bruno Pedroni3 , Kenneth Kreutz-Delgado2 , and Gert Cauwenberghs1,3
1

arXiv:1311.0966v3 [] 9 Dec 2013

2

Institute for Neural Computation, UCSD, La Jolla, CA, USA
Electrical and Computer Engineering Department, UCSD, La Jolla, CA, USA
3
Department of Bioengineering, UCSD, La Jolla, CA, USA
March 12, 2018

Abstract

posed of leaky I&F neurons with STDP synapses to learn a
generative model of the MNIST hand-written digit dataset,
and by testing it in recognition, generation and cue integration tasks.
Our results contribute to a machine learning-driven approach for synthesizing networks of spiking neurons capable of carrying out practical, high-level functionality.

Restricted Boltzmann Machines (RBMs) and Deep Belief
Networks have been demonstrated to perform efficiently in
a variety of applications, such as dimensionality reduction,
feature learning, and classification. Their implementation
on neuromorphic hardware platforms emulating large-scale
networks of spiking neurons can have significant advantages from the perspectives of scalability, power dissipation and real-time interfacing with the environment. However the traditional RBM architecture and the commonly
used training algorithm known as Contrastive Divergence
(CD) are based on discrete updates and exact arithmetics
which do not directly map onto a dynamical neural substrate. Here, we present an event-driven variation of CD to
train a RBM constructed with Integrate & Fire (I&F) neurons, that is constrained by the limitations of existing and
near future neuromorphic hardware platforms. Our strategy
is based on neural sampling, which allows us to synthesize
a spiking neural network that samples from a target Boltzmann distribution. The recurrent activity of the network replaces the discrete steps of the CD algorithm, while Spike
Time Dependent Plasticity (STDP) carries out the weight
updates in an online, asynchronous fashion.
We demonstrate our approach by training an RBM com∗

1

Introduction

Machine learning algorithms based on stochastic neural
network models such as RBMs and deep networks are currently the state-of-the-art in several practical tasks [28, 4].
The training of these models requires significant computational resources, and is often carried out using powerhungry hardware such as large clusters [35] or graphics processing units [5]. Their implementation in dedicated hardware platforms can therefore be very appealing from the
perspectives of power dissipation and of scalability.
Neuromorphic Very Large Scale Integration (VLSI) systems exploit the physics of the device to emulate very
densely the performance of biological neurons in a realtime fashion, while dissipating very low power [39, 30].
The distributed structure of RBMs suggests that neuromorphic VLSI circuits and systems can become ideal candidates for such a platform. Furthermore, the communication

Electronic address: nemre@ucsd.edu; Corresponding author

1

in software and in hardware platforms.
In standard CD, weight updates are computed on the
basis of alternating, feed-forward propagation of activities
[29]. In a neuromorphic implementation, this translates to
reprogramming the network connections and resetting its
state variables at every step of the training. As a consequence, it requires two distinct dynamical systems: one for
normal operation (i.e. testing), the other for training, which
is highly impractical. To overcome this problem, we train
the neural RBMs using an online adaptation of CD. We
exploit the recurrent structure of the network to mimic the
discrete “construction” and “reconstruction” steps of CD in
a spike-driven fashion, and Spike Time Dependent Plasticity (STDP) to carry out the weight updates. Each sample
(spike) of each random variable (neuron) causes synaptic
weights to be updated. We show that, over longer periods of
time, these microscopic updates behave like a macroscopic
CD weight update. Compared to standard CD, no additional
programming overhead is required during the training steps,
and both testing and training take place in the same dynamical system.
Because RBMs are generative models, they can act simultaneously as classifiers, content-addressable memories,
and carry out probabilistic inference. We demonstrate these
features in a MNIST hand-written digit task [37], using an
RBM network consisting of 824 “visible“ neurons and 500
“hidden” neurons. The spiking neural network was able
to learn a generative model capable of recognition performances with accuracies up to 91.9%, which is close to the
performance obtained using standard CD and Gibbs sampling, 93.6%.

between neuromorphic components is often mediated using asynchronous address-events [14] enabling them to be
interfaced with event-based sensors [38, 44, 42] for embedded applications, and to be implemented in a very scalable
fashion [31, 52, 50].
Currently, RBMs and the algorithms used to train them
are designed to operate efficiently on digital processors, using batch, discrete-time, iterative updates based on exact
arithmetic calculations. However, unlike digital processors,
neuromorphic systems compute through the continuoustime dynamics of their components, which are typically Integrate & Fire (I&F) neurons [30], rendering the transfer of
such algorithms on such platforms a non-trivial task. We
propose here a method to construct RBMs using I&F neuron models and to train them using an online, event-driven
adaptation of the Contrastive Divergence (CD) algorithm.
We take inspiration from computational neuroscience to
identify an efficient neural mechanism for sampling from
the underlying probability distribution of the RBM. Neuroscientists argue that brains deal with uncertainty in their
environments by encoding and combining probabilities optimally [18], and that such computations are at the core of
cognitive function [25]. While many mechanistic theories
of how the brain might achieve this exist, a recent neural
sampling theory postulates that the spiking activity of the
neurons encodes samples of an underlying probability distribution [20]. The advantage for a neural substrate in using such a strategy over the alternative one, in which neurons encode probabilities, is that it requires exponentially
fewer neurons. Furthermore, abstract model neurons consistent with the behavior of biological neurons can implement Markov Chain Monte Carlo (MCMC) sampling [7],
and RBMs sampled in this way can be efficiently trained using CD, with almost no loss in performance [45]. We identify the conditions under which a dynamical system consisting of I&F neurons performs neural sampling. These
conditions are compatible with neuromorphic implementations of I&F neurons, suggesting that they can achieve similar performance. The calibration procedure necessary for
configuring the parameters of the spiking neural network is
based on firing rate measurements, and so is easy to realize

2
2.1

Materials and Methods
Neural Sampling with Noisy I&F Neurons

We describe here the conditions under which a dynamical
system composed of I&F neurons can perform neural sampling. In has been proven that abstract neuron models consistent with the behavior of biological spiking neurons can
perform MCMC sampling of a Boltzmann distribution [7].
2

Two conditions are sufficient for this. First, the instantaneous firing rate of the neuron verifies:
(
0
if t − t0 < τr
,
(1)
ρ(u(t), t − t0 ) =
r(u(t)) t − t0 ≥ τr

synapses producing such PSPs is very difficult to realize in
hardware, when compared to first-order linear filters that result in “alpha”-shaped PSPs [17, 3]. This is because, in the
latter model, the synaptic dynamics are linear, such that a
single hardware synapse can be used to generate the same
current that would be generated by an arbitrary number of
synapses (see also next section). As a consequence, we will
use alpha-shaped PSPs instead of rectangular PSPs in our
models. The use of the alpha PSP over the rectangular PSP
is the major source of degradation in sampling performance,
as we will discuss in Sec. 2.2.

with r(u(t)) proportional to exp(u(t)), where u(t) is the
membrane potential and τr is an absolute refractory period
during which the neuron cannot fire. ρ(u(t), t − t0 ) describes the neuron’s instantaneous firing rate as a function
of u(t) at time t, given that the last spike occurred at t0 . The
average firing rate of this neuron model for stationary u(t)
is the sigmoid function:
ρ(u) = (τr + exp(−u))−1 .

Stochastic I&F Neurons. A neuron whose instantaneous
firing rate is consistent with Eq. (1) can perform neural
sampling. Eq. (1) is a generalization of the Poisson process to the case when the firing probability depends on
the time of the last spike (i.e. it is a renewal process),
and so can be verified only if the neuron fires stochastically [11]. Stochasticity in I&F neurons can be obtained
through several mechanisms, such as a noisy reset potential, noisy firing threshold, or noise injection [47]. The
first two mechanisms necessitate stochasticity in the neuron’s parameters, and therefore may require specialized circuitry. But noise injection in the form of background Poisson spike trains requires only synapse circuits, which are
present in many neuromorphic VLSI implementation of
spiking neurons [30, 3]. Furthermore, Poisson spike trains
can be generated self-consistently in balanced excitatoryinhibitory networks [55], or using finite-size effects and
neural mismatch [1].
We show that the abstract neuron model in Eq. (1) can be
realized in a simple dynamical system consisting of leaky
I&F neurons with noisy currents. The neuron’s membrane
potential below firing threshold θ is governed by the following differential equation:

(2)

Second, the membrane potential of neuron i is equal to the
linear sum of its inputs:
ui (t) = bi +

N
X

wij zj (t), ∀i = 1, ..., N,

(3)

j=1

where bi is a constant bias, and zj (t) represents the presynaptic spike train produced by neuron j and is set to 1
for a duration τr after the neuron has spiked. The terms
wij zj (t) are identified with the time course of the Post–
Synaptic Potential (PSP), i.e. the response of the membrane
potential to a pre-synaptic spike. The two conditions above
define a neuron model, to which we refer as the “abstract
neuron model”. The network then samples from a Boltzmann distribution:

1
exp − E(z1 , ..., zk ) , with
Z
X
1X
E(z1 , ..., zk ) = −
Wij zi zj −
bi zi ,
2
p(z1 , ...., zk ) =

ij

(4)

i

where Z is the partition function, and E(z1 , ..., zk ) can be
interpreted as an energy function [26].
An important fact of the abstract neuron model is that, according to the dynamics of zj (t), the PSPs are “rectangular”
and non-additive. The implementation of a large number of

C

d
ui = −gL ui + Ii (t) + σξ(t),
dt

ui (t) ∈ (−∞, θ),
(5)

where C is a membrane capacitance, ui is the membrane
potential of neuron i, gL is a leak conductance, σξ(t) is a
3

ρ(u0 )[Hz]

white noise term of amplitude σ (which can for example be
generated by background activity), Ii (t) its synaptic current
and θ is the neuron’s firing threshold. When the membrane
potential reaches θ, an action potential is elicited. After a
spike is generated, the membrane potential is clamped to
the reset potential urst for a refractory period τr .
In the case of the neural RBM, the currents Ii (t) depend
on the layer the neuron is situated in. For a neuron i in layer
v
N

125
3

0

0

4

u0

0
120

(6)

where Iid (t) is a current representing the data (i.e. the external input), I v is the feedback from the hidden layer activity
and the bias, and the q’s are the respective synaptic weights.
For a neuron j in layer h,

60

u0

0

30

θ−urst
σ =0.1

u0

0

Figure 1: Transfer curve of a leaky I&F neuron for three
different parameter sets where u0 = gIL , and τ1r = 250[Hz]
(dashed grey). In this plot, σV is varied to produce different
rst
ratios θ−u
σV . The three plots above shows that the fit with
the sigmoid function (solid black) improves as the ratio decreases.

Ij (t) =Ijh (t),
N

250

125
0

θ−urst
σ =0.3

125

j=1

v
X
d
qvij vi (t) + qbj bhj (t),
τsyn Ijh = − Ijh +
dt

250

250

Ii (t) =Iid (t) + Iiv (t),
h
X
d
τsyn Iiv = − Iiv +
qhji hj (t) + qbi bvi (t),
dt

θ−urst
σ =5.0

(7)

i=1

where I h is the feedback from the visible layer, and bv (t)
and bh (t) are Poisson spike trains implementing the bias.
The dynamics of I h and I v correspond to a first-order linear
filter, so each incoming spike results in PSPs that rise and
decay exponentially (i.e. alpha-PSP) [23].
Can this neuron neuron verify the conditions required for
neural sampling? The membrane potential is already assumed to be equal to the sum of the PSPs as required by
neural sampling. So to answer the above question we only
need to verify whether Eq. (1) holds. Eq. (5) is a Langevin
equation which can be analyzed using the Fokker-Planck
equation [22]. The solution to this equation provides the
neuron’s input/output response, i.e. its transfer curve (for a
review, see [48]):
!−1
Z θ−u0
σV
√
dx exp(x2 )(1 + erf(x))
,
ρ(u0 ) = τrf + τm π
urst −u0
σV

where erf is the error function (the integral of the normal
distribution), u0 = gI is the stationary value of the membrane potential when injected by a constant current I, urst
is the reset voltage, and σV2 (u) = σ 2 /(gL C).
According to Eq. (2), the condition for neural sampling
requires that the average firing rate of the neuron to be the
sigmoid function. Although the transfer curve of the noisy
I&F neuron Eq. (8) is not equal to the sigmoid function,
it was previously shown that with an appropriate choice of
parameters, the shape of this curve can be very similar to it
[40]. We observe that, for a given refractory period τr , the
rst
smaller ratio θ−u
in Eq. (5), the better the transfer curve
σV
rst
resembles a sigmoid function (Fig. 1) . With a small θ−u
σV ,
the transfer function of a neuron can be fitted to


1
exp(−Iβ) −1
ν(I) =
1+
,
(9)
τr
γτr
where β and γ are the parameters to be fitted. The choice

(8)
4

Neural mismatch can cause β and γ to differ from neuron
to neuron. From Eq. (9) and the linearity of the postsynaptic
currents I(t) in the weights, it is clear that this type of mismatch can be compensated by scaling the synaptic weights
and biases accordingly. The calibration of the parameters
γ and β quantitatively relate the spiking neural network’s
parameters to the RBM. In practice, this calibration step is
only necessary for mapping pre-trained parameters of the
RBM onto the spiking neural network.
Although we estimated the parameters of software simulated I&F neurons, parameter estimation based on firing
rate measurements were shown to be an accurate and reliable method for VLSI I&F neurons as well [43].

of the neuron model described in Eq. (5) is not critical for
neural sampling: A relationship that is qualitatively similar
to Eq. (8) holds for neurons with a rigid (reflective) lower
boundary [21] which is common in VLSI neurons, and for
I&F neurons with conductance-based synapses [46].
This result also shows that synaptic weights qvi , qhj ,
which have the units of charge are related to the RBM
weights Wij by a factor β −1 . To relate the neural activity to
the Boltzmann distribution, Eq. (4), each neuron is associated to a binary random variable which is assumed to take
the value 1 for a duration τr after the neuron has spiked,
and zero otherwise, similarly to [7]. The relation between
the random vector and the I&F neurons’ spiking activity is
illustrated in Fig. 3.

2.2

Calibration Protocol. In order to transfer the parameters
from the probability distribution Eq. (4) to those of the I&F
neurons, the parameters γ, β in Eq. (9) need to be fitted. An
estimate of a neuron’s transfer function can be obtained by
computing its spike rate when injected with different values
of constant inputs I. The refractory period τr is the inverse
of the maximum firing rate of the neuron, so it can be easily
measured by measuring the spike rate for very high input
current I. Once τr is known, the parameter estimation can
be cast into a simple linear regression problem by fitting
log(ρ(i)−1 −τr ) with βI +log(γ). Fig. 2 shows the transfer
curve when τr = 0 ms, which is approximately exponential
in agreement with Eq. (1).
The shape of the transfer curse is strongly dependent on
the noise amplitude. In the absence of noise, the transfer
curve is a sharp threshold function, which softens as the
amplitude of the noise is increased (Fig. 1). As a result,
both parameters γ and β are dependent on the variance of
the input currents from other neurons I(t). Since βq = w,
the effect of the fluctuations on the network is similar to
scaling the synaptic weights and the biases which can be
problematic. However, by selecting a large enough noise
amplitude σ and a slow enough input synapse time constant,
the fluctuations due to the background input are much larger
than the fluctuations due to the inputs. In this case, β and γ
remain approximately constant during the sampling.

Validation of Neural Sampling using I&F
neurons

The I&F neuron verifies Eq. (1) only approximately, and
the PSP model is different than the one of Eq. (3). Therefore, the following two important questions naturally arise:
how accurately does the I&F neuron-based sampler outlined above sample from a target Boltzmann distribution?
How well does it perform in comparison to an exact sampler, such as the Gibbs sampler? To answer these questions
we sample from several neural RBM consisting of 5 visible and 5 hidden units for randomly drawn weight and bias
parameters. At these small dimensions, the probabilities associated to all possible values of the random vector z can be
computed exactly. These probabilities are then compared to
those obtained through the histogram constructed with the
sampled events. To construct this histogram, each spike was
extended to form a box of length τr (as illustrated in Fig. 3),
the spiking activity was sampled at 1 kHz, and the occurrences of all the possible 210 states of the random vector z
were counted. We added 1 to the number of occurrences of
each state to avoid zero probabilities.
A common measure of similarity between two distributions is the Kullback-Leibler (KL) divergence:
D(p||q) =

X
i

5

pi log

pi
.
qi

νbias
σ

Mean firing rate of bias Poisson spike train
Noise amplitude

β
γ
τr
τsyn
τbr
gL
urst
C
θ
W
bv , bh
Nv , Nh

Exponential factor (fit)
Baseline firing rate (fit)
Refractory period
Time constant of recurrent, and bias synapses.
“Burn-in” time of the neural sampling
Leak conductance
Reset Potential
Membrane capacitance
Firing threshold
RBM weight matrix (∈ RNv ×Nh )
RBM bias for layer v and h
Number of visible and hidden units
in the RBM

Nc
2T

Number of class label units
Epoch duration

Tsim

Simulation time

τST DP
η

Learning time window
Learning rate

Table 1

6

all figures
all figures, except Fig. 1
Fig. 1 (left)
Fig. 1 (right)
Fig. 1 (bottom)
all figures
all figures
all figures
all figures
all figures
all figures
all figures
all figures
all figures
Fig. 4
Fig. 4
Fig. 4
Fig. 8,7
Fig. 9
Fig. 8,7,9
Fig. 4, 8,7
Fig. 9
Fig. 2
Fig. 4
Fig. 8
Fig. 9
Fig. 7 (testing)
Fig. 7 (learning)
Fig. 8
standard CD
event-driven CD

1000 Hz
3 · 10−11 A
2 · 10−11 A
3 · 10−10 A
1 · 10−9 A
2.044 · 109 A−1
8808 Hz
4 ms
4 ms
10 ms
1 nS
0V
10−12 F
100 mV
N (−.75, 1.5)
N (−1.5, .5)
5, 5
824, 500
834, 500
40
100 ms
300 ms
5s
1000 s
.2 s
.85 s
1.0 s
2000 s
4 ms
.1 · 10−2
3.2 · 10−2

If the distributions p and q are identical then D(p||q) = 0,
otherwise D(p||q) > 0.
The average KL divergence for 48 randomly drawn distributions after 1000 s of sampling time was 0.058 . This
result is not significantly different if the abstract neuron
model Eq. (1) with alpha PSPs is used, and in both cases
the KL divergence did not tend to zero as the number of
samples increased. The only difference in the latter neuron
model compared to the abstract neuron model of [7], which
tends to zero when sampling time tends to infinity, is the
PSP model. This indicates that the discrepancy is largely
due to the use of alpha-PSPs, rather than the approximation
of Eq. (1) with I&F neurons.
The standard sampling procedure used in RBMs is Gibbs
Sampling: the neurons in the visible layer are sampled simultaneously given the activities of the hidden neurons,
then the hidden neurons are sampled given the activities of
the visible neurons. This procedure is repeated a large number of times. For comparison with the neural sampler, the
duration of one Gibbs sampling iteration is identified with
one refractory period τr = 4 ms. At this scale, we observe
that the speed of convergence of the neural sampler is similar to that of the Gibbs sampler up to 104 ms, after which
the neural sampler plateaus above the D(p||q) = 10−2 line.
Despite the approximations in the neuron model and the
synapse model, these results show that in RBMs of this size,
the neural sampler consisting of I&F neurons sample from
a distribution that has the same KL divergence as the distribution obtained after 104 iterations of Gibbs sampling,
which is more than the typical number of iterations used for
MNIST hand-written digit tasks in the literature [27].

Figure 2: Transfer function of I&F neurons driven by background white noise (Eq. (5)). We measure the firing rate
of the neuron as a function of a constant current injection
to estimate ρ(u0 ), where for constant Iinj , u0 = Iinj /gL .
(Top) The transfer function of noisy I&F neurons in the absence of refractory period (ρ(u) = r(u), circles). We observe that ρ is approximately exponential over a wide range
of inputs, and therefore compatible with neural sampling.
Crosses show the transfer curve of neurons implementing
the abstract neuron Eq. (1), exactly. (Bottom) With an absolute refractory period the transfer function approximates
the sigmoid function. The firing rate saturates at 250 Hz
due to the refractory period chosen for the neuron. .

2.3

Neural Architecture for Learning a Model of
MNIST Hand-Written Digits

We test the performance of the neural RBM in a digit recognition task. We use the MNIST database, whose data samples consist of centered, gray-scale, 28×28-pixel images of
hand-written digits 0 to 9 [37]. The neural RBM’s network
architecture consisted of 2 layers, as illustrated in Fig. 5.
The visible layer was partitioned into 784 sensory neurons
7

Figure 4: (Left) Example probability distribution obtained by neural sampling of the RBM of Fig. 3. The bars are marginal
probabilities computed by counting the events [00000], [00001], · · · , [11110], [11111], respectively. PN S is the distribution
obtained by neural sampling and P is the exact probability distribution computed with Eq. (4). (Right) The degree to which
the sampled distribution resembles the target distribution is quantified by the KL divergence measured across 48 different
distributions, and the shadings correspond to its standard deviation. This plot also shows the KL divergence of the target
distribution sampled by Gibbs Sampling (PGibbs ), which is the common choice for RBMs. For comparison with the neural
sampler, we identified the duration of one Gibbs sampling iteration with one refractory period τr = 4 ms. The plot shows
that up to 104 ms, the two methods are comparable. After this, the KL divergence of the neural sampler tends to a plateau
due to the fact that neural sampling with our I&F neural network is approximate. In both figures, PN S,Abstract refers to the
marginal probability distribution obtained by using the abstract neuron model Eq. (1). In this case, the KL divergence is not
significantly different from the one obtained with the I&F neuron model-based sampler.

8

Hidden
(500 neurons)

Visible
(784 neurons)

Class
(40 neurons)

Data Layer
784+40 neurons

Figure 5: The RBM network consists of a visible and a
hidden layer. The visible layer is partitioned into 784 sensory neurons (vd ) and 40 class label neurons (vc ) for supervised learning. During data presentation, the activities
in the visible layer are driven by a data layer d, consisting
of a digit and its label (1 neuron per label). In the RBM,
the weight matrix between the visible layer and the hidden
layer is symmetric.
Figure 3: Neural Sampling in an RBM consisting of 10
stochastic I&F neurons, with 5 neurons in each layer. Each
neuron is associated to a binary random variable which take
values 1 during a refractory period τr after the neuron has
spiked (gray shadings). The variables are sampled at 1 kHz
to produce binary vectors that correspond to samples of the
joint distribution p(z). In this figure, only the membrane
potential and the samples produced by the first 5 neurons
are shown. The vectors inside the brackets are example
samples of the marginalized distribution p(z1 , z2 , z3 , z4 , z5 )
produced at the time indicated by the vertical lines. In the
RBM, there are no recurrent connections within a layer.

(vd ) and 40 class label neurons (vc ) for supervised learning. The pixel values of the digits were discretized to 2
values, with low intensity pixel values (p <= .5) mapped
to 10−5 and high intensity values (p > .5) mapped to 0.98.
A neuron i in d stimulated each neuron i in layer v, with
synaptic currents fi such that P (vi = 1) = ν(fi )τr = pi ,
where 0 ≤ pi ≤ 1 is the value of pixel i. The value fi
is calculated by invertingthe transfer
 function of the neus
−1
ron: fi = ν (s) = log γ−sγτr β −1 . Using this RBM,
classification is performed by choosing the most likely label given the input, under the learned model. This equals to
choosing the population of class neurons associated to the
same label that has the highest population firing rate.
To reconstruct a digit from a class label, the class neurons
belonging to a given digit are clamped to a high firing rate.
For testing the discrimination performance of an energybased model such as the RBM, it is common to compute
9

the free-energy F (vc ) of the class units [26], defined as:
X
exp(−F (vc )) =
exp(−E(vd , vc , h)),
(10)

and reconstruction phases. In practice, when the data set is
very large, weight updates are calculated using a subset of
data samples, or “mini-batches”. The above rule can then
be interpreted as a stochastic gradient descent [49]. Although the convergence properties of the CD rule are the
subject of continuing investigation, extensive software simulations show that the rule often converges to very good
solutions [29].
The main result of this paper is an online variation of the
CD rule for implementation in neuromorphic hardware. By
virtue of neural sampling the spikes generated from the visible and hidden units can be used to compute the statistics of
the probability distributions online (further details on neural
sampling in the Materials and Methods Sec. 2.1). Therefore
a possible neural mechanism for implementing CD is to use
synapses whose weights are governed by synaptic plasticity. Because the spikes cause the weight to update in an
online, and asynchronous fashion, we refer to this rule as
event-driven CD.
The weight update in event-driven CD is a modulated,
pair-based STDP rule:

vd ,h

and selecting vc such that the free-energy is minimized.
The spiking neural network is simulated using the BRIAN
simulator [24]. All the parameters used in the simulations
are provided in Tab. 1.

3
3.1

Results
Event-Driven Contrastive Divergence

A Restricted Boltzmann Machine (RBM) is a stochastic
neural network consisting of two symmetrically interconnected layers composed of neuron-like units - a set of visible units v and a set of hidden units h, but has no connections within a layer.
The training of RBMs commonly proceeds in two
phases. At first the states of the visible units are clamped
to a given vector from the training set, then the states of
the hidden units are sampled. In a second “reconstruction”
phase, the network is allowed to run freely. Using the statistics collected during sampling, the weights are updated in
a way that they maximize the likelihood of the data [29].
Collecting equilibrium statistics over the data distribution
in the reconstruction phase is often computationally prohibitive. The CD algorithm has been proposed to mitigate
this [29, 28]: the reconstruction of the visible units activity
is achieved by sampling them conditioned on the values of
the hidden units (Fig. 6). This procedure can be repeated k
times (the rule is then called CDk ), but relatively good convergence is obtained for the equilibrium distribution even
for one iteration. The CD learning rule is summarized as
follows:
∆wij = (hvi hj idata − hvi hj irecon ),

d
qij = g(t) STDPij (vi (t), hj (t))
dt

(12)

where g(t) ∈ R is a zero-mean global gating signal controlling the data vs. reconstruction phase, qij is the weight
of the synapse and vi (t) and hj (t) refer to the spike trains
of neurons vi and hj , respectively, which are represented
by a sum of Dirac delta pulses
P centered on the respective
hj (t) =
k∈Spi δ(t − tk ),
P spike times: vi (t) =
δ(t−t
)
where
Sp
and
Sp
are
the
set
of
the
spike
i
j
k
k∈Spj
times of the visible neuron i and hidden neuron j, respectively and δ(t) = 1 if t = 0 and 0 otherwise.
As opposed to the standard CD rule, weights are updated
after every occurrence of a pre-synaptic and post-synaptic
event. While this online approach slightly differentiates
it from standard CD, it is integral to a spiking neuromorphic framework where the data samples and weight updates cannot be stored. The weight update is governed by a
symmetric STDP rule with a symmetric temporal window

(11)

where vi and hj are the activities in the visible and hidden layers, respectively. This rule can be interpreted as a
difference of Hebbian and anti-Hebbian learning rules between the visible and hidden neurons sampled in the data
10

(b)

(a)

q
(data)

...
tpost-tpre

(reconstruction)

q

Figure 6: The standard Contrastive Divergence (CD)k procedure, compared to event-driven CD. (a) In standard CD, learning
proceeds iteratively by sampling in “construction” and “reconstruction” phases [29], which is impractical in a continuoustime dynamical system. (b) We propose a spiking neural sampling architecture that folds these updates on a continuous time
dimension through the recurrent activity of the network. The synaptic weight update follows a STDP rule modulated by a
zero mean signal g(t). This signal switches the behavior of the synapse from Long-Term Potentiation (LTP) to Long-Term
Depression (LTD), and partitions the training into two phases analogous to those of the original CD rule. The spikes cause
microscopic weight modifications, which on average behave as the macroscopic CD weight update. For this reason, the
learning rule is referred to as event-driven CD.

K(t) = K(−t), ∀t:
STDPij (vi (t), hj (t)) =vi (t)Ahj (t) + hj (t)Avi (t),
Z t
dsK(s − t)hj (s),
Ahj (t) =A
−∞
t

periods of the neurons [7]. In this work, we used the following modulation function g(t):

(13)

Z
Avi (t) =A

dsK(s − t)vi (s),


if mod(t, 2T ) ∈ (τbr , T )

1
g(t) = −1 if mod(t, 2T ) ∈ (T + τbr , 2T ) ,


0
otherwise

−∞

with A > 0 defining the magnitude of the weight updates.
In our implementation, updates are additive and weights can
change polarity.
Pairwise STDP with a global modulatory signal approximates CD. The modulatory signal g(t) switches the behavior of the synapse from LTP to LTD (i.e. Hebbian to
Anti-Hebbian). The temporal average of g(t) must vanish
to balance LTP and LTD, and must vary on much slower
time scales than the typical times scale of the network dynamics, denoted τbr , so that the network samples from its
stationary distribution when the weights are updated. The
time constant τbr corresponds to a “burn-in” time of MCMC
sampling and depends on the overall network dynamics and
cannot be computed in the general case. However, it is reasonable to assume τbr to be in the order of a few refractory

(14)

where mod is the modulo function and T is a time interval. The data is presented during the time intervals
(2iT, (2i + 1)T ), where i is a positive integer. With the
g(t) defined above, no weight update is undertaken during
a fixed period of time τbr . This allows us to neglect the transients after the stimulus is turned on and off (respectively in
the beginning of the data and reconstruction phases). In this
case and under further assumptions discussed below, the
event-driven CD rule can be directly compared with standard CD as we now demonstrate. The average weight up11

R0
br
with η := 2A T −τ
2T
−∞ d∆K(∆). Similar arguments apply to the averages in the time interval tr :

date during (0, 2T ) is:

d
qij i(0,2T ) =Cij + Rij ,
Z 0
dt
d∆K(∆)hvi (t)hj (t + ∆)itr = ηv̄i− h̄−
R
=
2A
ij
j .
T − τbr
−∞
(hvi (t)Ahj (t)itd + hhj (t)Avi (t)itd )
Cij =
2T
− −
T − τbr
Rij = −
(hvi (t)Ahj (t)itr + hhj (t)Avi (t)itr ), with v̄i h̄j := hvi (t)itr hhj (t + ∆)itr . The average update
2T
in (0, 2T ) then becomes:
(15)


d
+ +
− −
q
i
=
η
v̄
h
h̄
−
v̄
h̄
.
(18)
ij
(0,2T
)
where td = (τbr , T ) and tr = (T + τbr , 2T ) denote the
i j
i j
dt
intervals during the positive and negative phases of g(t),
Rb
1
According to Eq. (17), any symmetric temporal window
and h·i(a,b) = b−a
a dt·.
that is much shorter than T can be used. For simplicWe write the first average in Cij as follows:
ity, we choose an exponential temporal window K(∆) =
Z T Z t
exp(∆/τST DP ) with decay rate τST DP  T (Fig. 6b). In
1
dsK(s − t)vi (t)hj (s), this case, η = 2A T −τbr τST DP .
dt
hvi (t)Ahj (t)itd = A
2T
T − τbr τbr
−∞
The modulatory function g(t) partitions the training into
Z T Z 0
1
=A
dt
d∆K(∆)hvi (t)hj (t + ∆), several epochs of duration 2T . Each epoch consists of a
T − τbr τbr
−∞
LTP phase during which the data is presented (construcZ 0
tion), followed by a free-running LTD phase (reconstruc=A
d∆K(∆)hvi (t)hj (t + ∆)itd .
tion). The weights are updated asynchronously during the
−∞
time interval in which the neural sampling proceeds, and
(16)
Eq. (18) tells us that its average resembles Eq. (11). HowIf the spike times are uncorrelated the temporal averages
ever, it is different in two ways: the averages are taken over
become a product of the average firing rates of a pair of
one data and reconstruction phase rather than a mini-batch
visible and hidden neurons [23]:
of data samples and their reconstructions; and more importantly, the synaptic weights are updated during the data and
hvi (t)hj (t + ∆)itd = hvi (t)itd hhj (t + ∆)itd =: v̄i+ h̄+
j .
the reconstruction phase, whereas in the CD rule, updates
are carried out at the end of the reconstruction phase. In the
If we choose a temporal window that is much smaller than
derivation above the effect of the weight modification on the
T , and since the network activity is assumed to be stationnetwork during an epoch 2T was neglected for the sake of
ary in the interval (τbr , T ), we can write (up to a negligible
mathematical tractability. In the following, we verify that
error [32])
despite this approximation, the event-driven CD performs
nearly as well as standard CD in a commonly used benchZ 0
mark task.
hvi (t)Ahj (t)itd = Av̄i+ h̄+
d∆K(∆).
(17)
j
h

−∞

3.2
In the uncorrelated case, the second term in Cij contributes the same amount, leading to:
Cij =

Learning a generative model of hand-written
digits

We train the RBM to learn a generative model of the
MNIST handwritten digits using event-driven CD (see

ηv̄i+ h̄+
j .
12

Sec. 2.3 for details). For training, 20000 digits selected randomly from a training set consisting of 10000 digits were
presented in sequence, with an equal number of samples for
each digit.
The raster plots in Fig. 8 show the spiking activity of
each layer before and after learning for epochs of duration 100 ms. The top panel shows the population-averaged
weight. After training, the sum of the upwards and downward excursions of the average weight is much smaller than
before training, because the learning is near convergence.
The second panel shows the value of the modulatory signal
g(t). The third panel shows the input current (Id ) and the
current caused by the recurrent couplings (Ih ).
Two methods to estimate the overall classification performance of the neural RBM can be used. The first is by
neural sampling: the visible layer is clamped to the digit
only, and the network is run for 1s. The known label is then
compared with the positions of the group of class neurons
that had the highest population rate. The second method
is by minimizing free-energy: the neural RBMs parameters
are extracted, and for each data sample, the class neurons
with the lowest free-energy (See Materials and Methods) is
compared to the known label. In both cases, recognition
was tested for 1000 data samples that were not used during
the training. The results are summarized in Fig. 7.
As a reference we provide the best performance achieved
using the standard CD and one unit per class label (Nc =
10) (Fig. 7, table row 1), 93.6%. By mapping the learned
parameters to the neural RBM the recognition accuracy
reached 92.6%.
When training a neural RBM of I&F neurons using
event-driven CD, the recognition result was 91.9% (Fig. 7,
table row 2). The performance of this RBM obtained by
minimizing its free-energy was 90.8%. The learned parameters performed well for classification using the free-energy
calculation which suggests that the network learned a model
that is consistent with the mathematical description of the
RBM.
In an energy-based model like the RBM the free-energy
minimization should give the upper bound on the discrimination performance [26]. For this reason, the fact that

Neural Sampling Recognition Accuracy

Sampling Duration [s]

Standard CD
Event-driven CD
Event-driven CD (8 bits)
Event-driven CD (5 bits)

Accuracy
Neural Sampler
92.6%
91.9%
91.6%
89.4%

Accuracy
Free-energy
93.6%
90.8%
91.0%
89.2%

Figure 7: To test recognition accuracy, the RBMs are sampled using the I&F neuron-based sampler for up to 1 s. The
classification is read out by identifying the group of class
label neurons that had the highest activity. This experiment is run for RBM parameter sets obtained by standard
CD (black, CD) and event-driven CD (green, eCD). To test
the robustness of the RBM, it was run with parameters obtained by event-driven CD discretized to 8 and 5 bits. In
all scenarios, the accuracy after 50 ms of sampling was
above 80% and after 1 s the accuracies typically reached
their peak at around 91.9%. The dashed horizontal lines
show the recognition accuracy obtained by minimizing the
free-energy (see text). The fact that the eCD curve (solid
green) surpasses its free-energy minimization performance
suggests that the RBM learns a model that is tailored to the
I&F spiking neural network.

13

them, and infer a digit by combining partial evidence.
These features are clearly illustrated in the following experiment (Fig. 9). First the digit 3 is presented (i.e. layer vd is
driven by layer d) and the correct class label in vc activated.
Second, the neurons associated to class label 5 are clamped,
and the network generated its learned version of the digit.
Third, the right-half part of a digit 8 is presented, and the
class neurons are stimulated such that only 3 or 6 are able
to activate (the other class neurons are inhibited, indicated
by the gray shading). Because the stimulus is inconsistent
with 6, the network settled to 3 and reconstructed the left
part of the digit.
The latter part of the experiment illustrates the integration of information between several partially specified
cues, which is of interest for solving sensorimotor transformation or multi-modal sensory cue integration problems
[15, 18, 10]. This feature has been used for auditory-visual
sensory fusion in a spiking Deep Belief Network (DBN)
model [44]. There, the authors trained a DBN with visual
and auditory data, which learned to associate the two sensory modalities, very similarly to how class labels and visual data are associated in our architecture. Their network
was able to resolve a similar ambiguity as in our experiment
in Fig. 9, but using auditory inputs instead of a class label.
In the digit generation mode, the trained network had a
tendency to be globally bistable, whereby the layer vd completely deactivated layer h. Since all the interactions between vd and vc take place through the hidden layer, vc
could not reconstruct the digit. To avoid this, we added
two populations of I&F neurons that were wired to layers
v and h, respectively. The parameters of these neurons and
their couplings were tuned such that each layer was strongly
excited when it’s average firing rate fell below 5 Hz.

the recognition accuracy is higher when sampling as opposed to using the free-energy method may appear puzzling. However, this is possible because the neural RBM
does not exactly sample from the Boltzmann distribution,
as explained in Sec. 2.2. This suggests that event-driven
CD compensates for the discrepancy between the distribution sampled by the neural RBM and the Boltzmann distribution, by learning a model that is tailored to the spiking
neural network.
Excessively long training durations can be impractical
for real-time neuromorphic systems. Fortunately, the learning using event-driven CD is fast: Compared to the off-line
RBM training (250000 presentations, in mini-batches of
100 samples) the event-driven CD training succeeded with
a smaller number of data presentations (20000), which corresponded to 2000 s of simulated time. This suggests that
the training durations are achievable for real-time neuromorphic systems.
The choice of the number of class neurons Nc . Eventdriven CD underperformed in the case of 1 neuron per class
label (Nc = 10), which is the common choice for standard
CD and Gibbs sampling. This is because a single neuron firing at its maximum rate of 250 Hz cannot efficiently drive
the rest of the network without tending to induce spike-tospike correlations (e.g. synchrony), which is incompatible
with the assumptions made for sampling with I&F neurons
and event-driven CD. As a consequence, the generative
properties of the neural RBM degrade. This problem is
avoided by using several neurons per class label (in our case
four neurons per class label) because the synaptic weight
can be much lower to achieve the same effect, resulting in
smaller spike-to-spike correlations.

3.3

Neural parameters with finite precision. In hardware
systems, the parameters related to the weights and biases
cannot be set with floating-point precision, as can be done
in a digital computer. In current neuromorphic implementations the synaptic weights can be configured at precisions
of about 8 bits [56]. We characterize the impact of finiteprecision synaptic weights on performance by discretizing

Generative properties of the RBM

We test the neural RBM as a generative model of the
MNIST dataset of handwritten digits, using parameters obtained by running the event-driven CD.
In the context of the handwritten digit task, the RBM’s
generative property enables it to classify digits, generate
14

Figure 8: The spiking neural network learns a generative model of the MNIST dataset using the event-driven CD procedure.
(a) Learning curve, shown here up to 10000 samples. (b) Details of the training procedure, before and after training (20000
samples). During the first half of each .1 s epoch, the visible layer v is driven by the sensory layer. During this phase, the
gating variable g is 1, meaning that the synapses undergo LTP. During the second half of each epoch, the sensory stimulus is
removed, and g is set to −1, so the synapses undergo LTD. The top panels of both figures show the mean of the entries of the
weight matrix. The second panel shows the values of the modulatory signal g(t). The third panel shows the synaptic currents
of a visible neuron, where Ih is caused by the feedback from the hidden and the bias, and Id is the data. The timing of the
clamping (Id ) and g differ due to an interval τbr where no weight update is undertaken to avoid the transients (See Materials
and Methods). Before learning and during the reconstruction phase, the activity of the visible layer is random. But as learning
progresses, the activity in the visible layer reflects the presented data in the reconstruction phase. This is very well visible in
the layer class label neurons vc , whose activity persists after the sensory stimulus is removed. Although the firing rates of the
hidden layer neurons before training is high (average 113 Hz), this is only a reflection of the initial conditions for the recurrent
couplings W . In fact, at the end of the training, the firing rates in both layers becomes much sparser (average 9.31 Hz).

15

Figure 9: The recurrent structure of the network allows it to classify, reconstruct and infer from partial evidence. (a) Raster
plot of an experiment illustrating these features. Before time 0s, the neural RBM runs freely, with no input. Due to the
stochasticity in the network, the activity wanders from attractor to attractor. At time 0s, the digit 3 is presented (i.e. layer vd
is driven by d), activating the correct class label in vc ; At time t = .3 s, the class neurons associated to 5 are clamped to high
activity and the rest of the class label neurons are strongly inhibited, driving the network to reconstruct its version of the digit
in layer vd ; At time t = .6 s, the right-half part of a digit 8 is presented, and the class neurons are stimulated such that only 3
or 6 can activate (all others are strongly inhibited as indicated by the gray shading). Because the stimulus is inconsistent with
6, the network settles to a 3 and attempts to reconstruct it. The top figures shows the digits reconstructed in layer vd . (b) Digits
0 − 9, reconstructed in the same manner. Each row corresponds to a different, independent run. (c) Population firing rate
of the experiment presented in (a). During recognition, the network typically reaches equilibrium after about 10τr = 40 ms
(black bar).

16

the weight and bias parameters to 8 bits and 5 bits. The
set of possible weights were spaced uniformly in the interval (µ − 4.5σ, µ + 4.5σ), where µ, σ are the mean and the
standard deviation of the parameters across the network, respectively. The classification performance of MNIST digits degraded gracefully. In the 8 bit case, it degrades only
slightly to 91.6%, but in the case of 5 bits, it degrades more
substantially to 89.4%. In both cases, the RBM still retains
its discriminative power, which is encouraging for implementation in hardware neuromorphic systems.

4

chine learning tasks [4]. The ability to implement RBMs
with spiking neurons and train then using event-based CD
paves the way towards on-line training of DBNs of spiking
neurons [27].
We chose the MNIST handwritten digit task as a benchmark for testing our model. When the RBM was trained
with standard CD, it could recognize up to 926 out of 1000
of out-of-training samples. The MNIST handwritten digits
recognition task was previously shown in a digital neuromorphic chip [2], which performed at 89% accuracy, and
in a software simulated visual cortex model [19]. However, both implementations were configured using weights
trained off-line. A recent article showed the mapping of
off-line trained DBNs onto spiking neural network [44].
Their results demonstrated hand-written digit recognition
using neuromorphic event-based sensors as a source of input spikes. Their performance reached up to 94.1% using
leaky I&F neurons. The use of off-line CD combined with
an additional layer explains to a large extent their better
performance compared to ours. Our work extends [44] by
demonstrating an on-line training using synaptic plasticity,
testing its robustness to finite weight precision, and providing an interpretation of spiking activity in terms of neural
sampling.
To achieve the computations necessary for sampling
from the RBM, we have used the neural sampling framework [20], where each spike is interpreted as a sample of an
underlying probability distribution. [7] proved that abstract
neuron models consistent with the behavior of biological
spiking neurons can perform MCMC, and have applied it to
a basic learning task in a fully visible Boltzmann Machine.
We extended the neural sampling framework in three ways:
First, we identified the conditions under which a dynamical system consisting of I&F neurons can perform neural
sampling; Second, we verified that the sampling of RBMs
was robust to finite-precision parameters; Third, we demonstrated learning in a Boltzmann Machine with hidden units
using STDP synapses.
In the neural sampling framework, neurons behave
stochastically. This behavior can be achieved in I&F neurons using noisy input currents, created by a Poisson spike

Discussion

Neuromorphic systems are promising alternatives for largescale implementations of RBMs and deep networks, but the
common procedure used to train such networks, Contrastive
Divergence (CD), involves iterative, discrete-time updates
that do not straightforwardly map on a neural substrate. We
solve this problem in the context of the RBM with a spiking neural network model that uses the recurrent network
dynamics to compute these updates in a continuous-time
fashion. We argue that the recurrent activity coupled with
STDP dynamics implements an event-driven variant of CD.
Using event-driven CD, the network connectivity remains
unchanged during training and testing, enabling the system
to learn in an on-line fashion, while being able to carry out
functionally relevant tasks such as recognition, data generation and cue integration.
The CD algorithm can be used to learn the parameters of probability distributions other than the Boltzmann
distribution (even those without any symmetry assumptions). Our choice for the RBM, whose underlying probability distribution is a special case of the Boltzmann distribution, is motivated by the following facts: They are
universal approximators of discrete distributions [36]; the
conditions under which a spiking neural circuit can naturally perform MCMC sampling of a Boltzmann distribution
were previously studied [40, 7]; and RBMs form the building blocks of many deep learning models such as DBNs,
which achieve state-of-the-art performance in many ma17

each neuron. Sharing one synapse circuit per pair of neurons can solve this problem. This may be impractical due
to the very large number of synapse circuits in the network,
but may be less problematic when using Resistive RandomAccess Memorys (RRAMs) (also called memristors) crossbar arrays to emulate synapses [12, 34, 51].RRAM are a
new class of nanoscale devices whose current-voltage relationship depends on the history of other electrical quantities
[53], and so act like programmable resistors. Because they
can conduct currents in both directions, one RRAM circuit
can be shared between a pair of neurons. A second problem is the number of recurrent connections. Even our RBM
of modest dimensions involved almost 2 million synapses,
which is impractical in terms of bandwidth and weight storage. Even if a very high number of weights are zero, the
connections between each pair of neurons must exist in order for a synapse to learn such weights. One possible solution is to impose sparse connectivity between the layers
[54, 41]. This remains to be tested in our model.

train. Spike trains with Poisson-like statistics can be generated with no additional source of noise, for example by the
following mechanisms: balanced excitatory and inhibitory
connections [55], finite-size effects in a large network, and
neural mismatch [1]. The latter mechanism is particularly
appealing, because it benefits from fabrication mismatch
and operating noise inherent to neuromorphic implementations [9].
Other groups have also proposed to use I&F neuron models for computing the Boltzmann distribution. [40] have
shown that noisy I&F neurons’ activation function is approximately sigmoidal as required by the Boltzmann machine, and have devised a scheme whereby a global inhibitory rhythm drives the network to generate samples
of the Boltzmann distribution. [44] have demonstrated a
deep belief network of I&F neurons that was trained offline, using standard CD and tested it using the MNIST
database. Independently and simultaneously to this work,
[46] demonstrated that conductance-based I&F neurons in
a noisy environment are compatible with neural sampling
as described in [7]. Similarly, [46] find that the choice of
non-rectangular PSPs and the approximations made by the
I&F neurons are not critical to the performance of the neural sampler. Our work extends all of those above by providing an online, STDP-based learning rule to train RBMs
sampled using I&F neurons.

Outlook: A custom learning rule. Our method combines I&F neurons that perform neural sampling and the
CD rule. Although we showed that this leads to a functional model, we do not know whether event-driven CD is
optimal in any sense. This is partly due to the fact that CDk
is an approximate rule [29], and it is still not entirely understood why it performs so well, despite extensive work
in studying its convergence properties [8]. Furthermore,
the distribution sampled by the I&F neuron does not exactly correspond to the Boltzmann distribution, and the average weight updates in event-driven CD differ from those
of standard CD, because in the latter they are carried out at
the end of the reconstruction step.
A very attractive alternative is to derive a custom synaptic plasticity rule that minimizes some functionally relevant quantity (such as Kullback-Leibler divergence or Contrastive Divergence), given the encoding of the information in the I&F neuron [16, 6]. A similar idea was recently pursued in [6], where the authors derived a tripletbased synaptic learning rule that minimizes an upper bound
of the Kullback-Leibler divergence between the model and

Applicability to neuromorphic hardware. Neuromorphic systems are sensible to fabrication mismatch and operating noise. Fortunately, the mismatch in the synaptic
weights and the activation function parameters γ and β are
not an issue if the biases and the weights are learned, and
the functionality of the RBM is robust to small variations in
the weights caused by discretization. These two findings are
encouraging for neuromorphic implementations of RBMs.
However, at least two conceptual problems of the presented
RBM architecture must be solved in order to implement
such systems on a large-scale. First, the symmetry condition required by the RBM does not necessarily hold. In
a neuromorphic device, the symmetry condition is impossible to guarantee if the synapse weights are stored locally at
18

and Y. Bengio. Theano: a CPU and GPU math expression compiler. In Proceedings of the Python for
Scientific Computing Conference (SciPy), volume 4,
2010.

the data distributions. Interestingly, their rule had a similar
global signal that modulates the learning rule, as in eventdriven CD, although the nature of this resemblance remains
to be explored. Such custom learning rules can be very beneficial in guiding the design of on-chip plasticity in neuromorphic VLSI and RRAM nanotechnologies, and will be
the focus of future research.

[6] Johanni Brea, Walter Senn, and Jean-Pascal Pfister.
Matching recall and storage in sequence learning with
spiking neural networks. The Journal of Neuroscience, 33(23):9565–9575, 2013.

Acknowledgments
[7] L. Buesing, J. Bill, B. Nessler, and W. Maass. Neural
dynamics as sampling: A model for stochastic computation in recurrent networks of spiking neurons. PLoS
Computational Biology, 7(11):e1002211, 2011.

This work was partially funded by the National Science
Foundation (NSF EFRI-1137279), the Office of Naval Research (ONR MURI 14-13-1-0205), and the Swiss National
Science Foundation (PA00P2_142058). We thank all the
anonymous reviewers of a previous version of this article
for their constructive comments.

[8] Miguel A Carreira-Perpinan and Geoffrey E Hinton.
On contrastive divergence learning. In Artificial Intelligence and Statistics, volume 2005, page 17, 2005.

[1] D. Amit and N. Brunel. Model of global spontaneous
activity and local structured activity during delay periods in the cerebral cortex. Cerebral Cortex, 7:237–
252, 1997.

[9] E. Chicca and S. Fusi. Stochastic synaptic plasticity in deterministic aVLSI networks of spiking neurons. In Frank Rattay, editor, Proceedings of the World
Congress on Neuroinformatics, ARGESIM Reports,
pages 468–477, Vienna, 2001. ARGESIM/ASIM Verlag.

[2] J. Arthur, P. Merolla, F. Akopyan, R. Alvarez, A. Cassidy, S. Chandra, S. Esser, N. Imam, W. Risk, DBD
Rubin, et al. Building block of a programmable neuromorphic substrate: A digital neurosynaptic core.
In Neural Networks (IJCNN), The 2012 International
Joint Conference on, pages 1–8. IEEE, 2012.

[10] D. Corneil, D. Sonnleithner, E. Neftci, E. Chicca,
M. Cook, G. Indiveri, and R. Douglas. Function approximation with uncertainty propagation in a VLSI
spiking neural network. In International Joint Conference on Neural Networks, IJCNN 2012, pages 2990–
2996. IEEE, 2012.

[3] C. Bartolozzi and G. Indiveri. Synaptic dynamics
in analog VLSI. Neural Computation, 19(10):2581–
2603, Oct 2007.

[11] D.R. Cox. Renewal theory, volume 1. Methuen London, 1962.

References

[4] Y. Bengio. Learning deep architectures for ai. Foundations and Trends R in Machine Learning, 2(1):1–127,
2009.

[12] Jose M Cruz-Albrecht, Timothy Derosier, and
Narayan Srinivasa. A scalable neural chip with
synaptic electronics using cmos integrated memristors. Nanotechnology, 24(38):384011, 2013.

[5] J. Bergstra, O. Breuleux, F. Bastien, P. Lamblin,
R. Pascanu, G. Desjardins, J. Turian, D. Warde-Farley,

[13] George E Dahl, Tara N Sainath, and Geoffrey E Hinton. Improving deep neural networks for lvcsr using
19

[23] W. Gerstner and W. Kistler. Spiking Neuron Models.
Single Neurons, Populations, Plasticity. Cambridge
University Press, 2002.

rectified linear units and dropout. In Proc. ICASSP,
2013.
[14] S.R. Deiss, R.J. Douglas, and A.M. Whatley. A pulsecoded communications infrastructure for neuromorphic systems. In W. Maass and C.M. Bishop, editors, Pulsed Neural Networks, chapter 6, pages 157–
78. MIT Press, 1998.

[24] D. Goodman and R. Brette. Brian: a simulator for
spiking neural networks in Python. Frontiers in Neuroinformatics, 2, 2008.
[25] T.L. Griffiths, N. Chater, C. Kemp, A. Perfors, and
J. B. Tenenbaum. Probabilistic models of cognition:
exploring representations and inductive biases. Trends
in cognitive sciences, 14(8):357–364, 2010.

[15] S. Deneve, PE Latham, and A. Pouget. Efficient computation and cue integration with noisy population
codes. Nature Neuroscience, 4(8):826–831, 2001.

[26] Simon S Haykin. Neural networks: A comprehensive
foundation, 1999.

[16] Sophie Deneve. Bayesian spiking neurons I: Inference. Neural Computation, 20(1):91–117, 2008.

[27] G.E. Hinton, S. Osindero, and Y.W. Teh. A fast learning algorithm for deep belief nets. Neural computation, 18(7):1527–1554, 2006.

[17] A. Destexhe, Z.F. Mainen, and T.J. Sejnowski. Methods in Neuronal Modelling, from ions to networks,
chapter Kinetic Models of Synaptic Transmission,
pages 1–25. MIT Press, 1998.

[28] G.E. Hinton and R.R. Salakhutdinov. Reducing the
dimensionality of data with neural networks. Science,
313(5786):504–507, 2006.

[18] K. Doya, S. Ishii, A. Pouget, and R.P.N. Rao.
Bayesian Brain Probabilistic Approaches to Neural
Coding. MIT Press, 2007.

[29] Geoffrey E Hinton. Training products of experts by
minimizing contrastive divergence. Neural computation, 14(8):1771–1800, 2002.

[19] C. Eliasmith, T.C. Stewart, X. Choo, T. Bekolay,
T. DeWolf, Y. Tang, and D. Rasmussen. A largescale model of the functioning brain.
Science,
338(6111):1202–1205, 2012.

[30] G. Indiveri, B. Linares-Barranco, T.J. Hamilton,
A. van Schaik, R. Etienne-Cummings, T. Delbruck, S.-C. Liu, P. Dudek, P. Häfliger, S. Renaud,
J. Schemmel, G. Cauwenberghs, J. Arthur, K. Hynna,
F. Folowosele, S. Saighi, T. Serrano-Gotarredona,
J. Wijekoon, Y. Wang, and K. Boahen. Neuromorphic
silicon neuron circuits. Frontiers in Neuroscience,
5:1–23, 2011.

[20] J. Fiser, P. Berkes, G. Orbán, and M. Lengyel. Statistically optimal perception and learning: from behavior
to neural representations: Perceptual learning, motor
learning, and automaticity. Trends in cognitive sciences, 14(3):119, 2010.

[31] S. Joshi, S. Deiss, M. Arnold, J. Park, T. Yu, and
G. Cauwenberghs. Scalable event routing in hierarchical neural array architecture with global synaptic connectivity. In Cellular Nanoscale Networks and Their
Applications (CNNA), 2010 12th International Workshop on, pages 1–6. IEEE, 2010.

[21] S. Fusi and M. Mattia. Collective behavior of networks with linear (VLSI) integrate and fire neurons.
Neural Computation, 11:633–52, 1999.
[22] Crispin W Gardiner. Handbook of stochastic methods.
2012.
20

complete sparse learning.
19(9):2301–2352, 2007.

[32] R. Kempter, W. Gerstner, and J.L. Van Hemmen.
Intrinsic stabilization of output rates by spike-based
hebbian learning. Neural Computation, 13(12):2709–
2741, 2001.

Neural Computation,

[42] E. Neftci, J. Binas, U. Rutishauser, E. Chicca, G. Indiveri, and R. J. Douglas. Synthesizing cognition in neuromorphic electronic systems. Proceedings of the National Academy of Sciences, 110(37):E3468–E3476,
June 2013.

[33] David F Kerridge. Inaccuracy and inference. Journal
of the Royal Statistical Society. Series B (Methodological), pages 184–194, 1961.
[34] Duygu Kuzum, Rakesh GD Jeyasingh, Byoungil Lee,
and H-S Philip Wong. Nanoelectronic programmable
synapses based on phase change materials for braininspired computing. Nano letters, 12(5):2179–2186,
2011.

[43] E. Neftci, B. Toth, G. Indiveri, and H. Abarbanel. Dynamic state and parameter estimation applied to neuromorphic systems. Neural Computation,
24(7):1669–1694, July 2012.
[44] P O’Connor, D. Neil, S.-C. Liu, T. Delbruck, and
M. Pfeiffer. Real-time classification and sensor fusion with a spiking deep belief network. Frontiers in
Neuroscience, 7(178), 2013.

[35] Quoc V Le, Marc’Aurelio Ranzato, Rajat Monga,
Matthieu Devin, Kai Chen, Greg S Corrado, Jeff
Dean, and Andrew Y Ng. Building high-level features using large scale unsupervised learning. arXiv
preprint arXiv:1112.6209, 2011.

[45] B. Pedroni, S. Das, E. Neftci, K. Kreutz-Delgado,
and G. Cauwenberghs. Neuromorphic adaptations of
restricted boltzmann machines and deep belief networks. In International Joint Conference on Neural
Networks, IJCNN 2013, 2013. (accepted).

[36] Nicolas Le Roux and Yoshua Bengio.
Representational power of restricted boltzmann machines
and deep belief networks. Neural Computation,
20(6):1631–1649, 2008.

[46] Mihai A Petrovici, Johannes Bill, Ilja Bytschok, Johannes Schemmel, and Karlheinz Meier. Stochastic
inference with deterministic spiking neurons. arXiv
preprint arXiv:1311.3211, 2013.

[37] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.
Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324,
1998.

[47] Hans E Plesser and Wulfram Gerstner. Noise in
integrate-and-fire neurons: From stochastic input to
escape rates. Neural Computation, 12(2):367–384,
2000.

[38] S.-C. Liu and T. Delbruck. Neuromorphic sensory systems. Current Opinion in Neurobiology,
20(3):288–295, 2010.
[39] C.A. Mead. Analog VLSI and Neural Systems.
Addison-Wesley, Reading, MA, 1989.

[48] A. Renart, P. Song, and X.-J. Wang. Robust spatial working memory through homeostatic synaptic
scaling in heterogeneous cortical networks. Neuron,
38:473–485, May 2003.

[40] P. Merolla, T. Ursell, and J. Arthur. The thermodynamic temperature of a rhythmic spiking network.
ArXiv e-prints, September 2010.

[49] Herbert Robbins and Sutton Monro. A stochastic approximation method. The Annals of Mathematical
Statistics, pages 400–407, 1951.

[41] Joseph F Murray and Kenneth Kreutz-Delgado. Visual recognition and inference using dynamic over21

[50] J. Schemmel, D. Brüderle, A. Grübl, M. Hock,
K. Meier, and S. Millner. A wafer-scale neuromorphic hardware system for large-scale neural modeling.
In International Symposium on Circuits and Systems,
ISCAS 2010, pages 1947–1950. IEEE, 2010.
[51] Teresa Serrano-Gotarredona, Timothée Masquelier,
Themistoklis Prodromakis, Giacomo Indiveri, and
Bernabe Linares-Barranco. Stdp and stdp variations
with memristors for spiking neuromorphic learning
systems. Frontiers in neuroscience, 7, 2013.
[52] R. Silver, K. Boahen, S. Grillner, N. Kopell, and K.L.
Olsen. Neurotech for neuroscience: unifying concepts, organizing principles, and emerging tools. J.
Neurosci., 27(44):11807, 2007.
[53] Dmitri B Strukov, Gregory S Snider, Duncan R Stewart, and R Stanley Williams. The missing memristor
found. Nature, 453(7191):80–83, 2008.
[54] Yichuan Tang and Chris Eliasmith. Deep networks
for robust visual recognition. In Proceedings of the
27th International Conference on Machine Learning
(ICML-10), pages 1055–1062, 2010.
[55] C. van Vreeswijk and H. Sompolinsky. Chaos in
neuronal networks with balanced excitatory and inhibitory activity. Science, 274(5293):1724–1726, December 1996.
[56] T. Yu and G. Cauwenberghs. Log-domain timemultiplexed realization of dynamical conductancebased synapses. In International Symposium on Circuits and Systems, (ISCAS), 2010, pages 2558 –2561.
IEEE, June 2010.

22

