SparseProp: Efficient Event-Based Simulation and Training of Sparse Recurrent Spiking Neural Networks
Authors: Rainer Engelken
Venue: NeurIPS 2023 poster (2023)
Keywords: spiking networks
event-based simulation
sparse networks
backpropagation
algorithm
neuroscience
TL;DR: We propose a novel event-based algorithm to simulate and train spiking neural networks, reducing computational cost from N to log(N) per network spike for sparse spiking networks.

Abstract: Spiking Neural Networks (SNNs) are biologically-inspired models that are capable of processing information in streams of action potentials. However, simulating and training SNNs is computationally expensive due to the need to solve large systems of coupled differential equations. In this paper, we propose a novel event-based algorithm called SparseProp for simulating and training sparse SNNs. Our algorithm reduces the computational cost of both forward pass and backward pass operations from O(N) to O(log(N)) per network spike, enabling numerically exact simulations of large spiking networks and their efficient training using backpropagation through time. By exploiting the sparsity of the network, SparseProp avoids iterating through all neurons at every spike and uses efficient state updates. We demonstrate the effectiveness of SparseProp for several classical integrate-and-fire neuron models, including simulating a sparse SNN with one million LIF neurons, which is sped up by more than four orders of magnitude compared to previous implementations. Our work provides an efficient and exact solution for training large-scale spiking neural networks and opens up new possibilities for building more sophisticated brain-inspired models.

SparseProp: Efficient Event-Based Simulation and Training of Sparse Recurrent Spiking Neural Networks
Rainer Engelken
Zuckerman Mind Brain Behavior Institute, Columbia University, New York, USA Max Planck Institute for Dynamics and Self -rfa--ngaazzainn,, Göttingen, Germany Göttingen Campus Institute for Dynamics of Biological Networks, Göttingen, Germany re2365@columbaa.dd
Abstract
Spiking Neural Networks (SNNs) are biologically-inspree models that capable of processing information in streams action potentials. However, simulating and training SNNs is computationally expensive due to the need to solve large systems of coupled differential equations. In this paper, we introduce SparseProp. a novel event-based algorithm for simulating and training sparse SNNs. Our algorithm reduces the computational cost of both the forward and backward pass operations from O(N) O(log(N)) per network spike, thereby enabling numerically exact sim- ulations of large spiking networks and their efficient training using backpropagation through time. By leveraging the sparsity of the network, SparseProp eliminates the need to iterate through all neurons at each spike, employing efficient state updates instead. We demonstrate the efficacy of SparseProp across several classical integrate-and-finn neuron models, including a simulation of sparse SNN with one million LIF neurons. This results speed-up exceeding four orders of magnitude relative to previous event- based implementations. Our work provides an efficient and exact solution for training large-scale spiking neural networks and opens up new possibilities for building more sophisticated brain-inspired models.
1 Introduction
The cortex processes information via streams action potentials commonly referred to as spikes that propagate within and between layers of recurrent neural networks. Spiking neural networks (SNNs) provide a more biologically plausible description of neuronal activity, capturing in contrast rate networks membrane potentials and temporal sequences of spikes. Moreover, SNNs promise solutions to the high energy consumption and CO2 emissions of deep networks [ [], as well as the
ability transmit information through precise spike timing 3 5).
Unlike commonly used iterative ODE solvers that necessitate time discretization, event-based sim- ulations resolve neural network dynamics precisely between spike events. This approach ensures machine-level precision in simulations, thereby alleviating the necessity for verification of result robustness with respect to time step size At and mitigating numerical issues, for instance in the vicinity of synchronous regimes 6 However, major drawback of event-based simulations is the computational cost of iterating through all neurons every network spike time.
To tackle this challenge, we introduce SparseProp, an novel event- algorithm designed for simu- lating and training sparse SNNs. Our algorithm reduces the computational cost of both the forward and backward pass from O(N) to (llgg(N) per network spike. This efficiency enables numerically exact simulations of large spiking networks and their efficient training using backpropagation through
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

time. By exploiting network sparsity and utilizing a change of variable to represent neuron states as times the next spikes on binary heap, SparseProp avoids iterating through all neurons every spike and employs efficient state updates.
We demonstrate the utility of PparseProp by applying it to three popular spiking neuron models. While the current implementation is for univariate neuron models, we discuss potential extension to more detailed multivariate models later.
Our contributions include:
Introducing SparseProp. novel and efficient algorithm for numerically exact event-based simulations of recurrent SNNs (section 2 and appendix A for minimal example code).
Conducting a numerical and analytical scaling analysis that compares the computational cost of our proposed algorithm to conventional event-based spiking network simulations (section 3 Fig and table [
Providing concrete implementations of the algorithm for recurrent networks of leaky
integrate-and-fire neurons and quadratic integrate-and-fire neurons (section and here)
Extending the algorithm to neuron models that lack an analytical solution for the next spike time using Chebyshev polynomials (section 5 and appendix B
* Extending the algorithm to heterogeneous spiking networks (section 6
In summary, our proposed SparseProp algorithm offers promising approach for simulating and training SNNs with both machine precision and practical applicability. It enables the simulation and training of large-scale SNNs with significantly reduced computational costs, paving the way for the advancement of brain-inspirdd models.
2 Algorithm for Efficient Event- Based Simulations in Recurrent Networks
We consider the dynamics of a spiking recurrent neural network of N neurons that is described by a system of coupled differential equations [7 10
12
13
14)
dV.
= F(V) + 1 Σ Jig h(t 1.)))
(1)
m dt
Here, the rate of change of the membrane potential Vi depends on its internal dynamics F(V,), an external input Iext and the recurrent input Σ Jiy h(t - e. Tm is the membrane time constant.
J.8
Thus, when presynaptic neuron j spikes at time 164) the hth postsynaptic neuron is driven by by some temporal kernel h(T) and coupling strength Jtj.
For event-based simulations, the dynamics of Eq aare are evolved from network spike to network spike, usually based on an analytical solution of the membrane potential Vt(t) as function of time t, instead of using iterative ODE solvers like the Runge- Kutta methods. A key aspect of this method for spiking networks is that the spike times are not confined to the time steps of a numerical solver, allowing them to be obtained with machine precision. To simulate the entire network dynamics this way involves four simple steps:
1. Find the next spiking neuron j and spike time 1y-.
2. Evolve the state V42 of every neuron to the next network spike time Lj+.
3. Update the neurons postsynaptic to j* using the weight Jıj.
4. Reset the state of the spiking neuron V,, Vreset-
Such event-based simulation schemes are usually used for neuron models that have an analytical solution of the individual neurons' dynamics Vi(t) [10, [4] [5]
16 I7 [8]. The precision of the solution of event-based simulation then limited only by the machine precision of the computer. In contrast, solving spiking network dynamics with an iterative ODE solver requires asserting that results are not dependent on the integration step size AL of the integration scheme. The total accumulated error usually scale with O (()t)),, where p is the order of the ODE solver. However, both iterating

---

A
B
PTC g(
voltage
22
-1
0.00
0.25
0.50
0.75
1.00
0
phase +
D
46
phase e
reset
1/
4.0
4$
0.00
0.25
0.50
0.75
1.00
Peeset
Prseet
4mn
time (s)
Figure 1: Change of reference frame allows efficient evolution of network state. A We map all membrane potentials (1) to phases (1) that evolve linearly between spikes with phase velocity w [10] 15) B At spike times t.. the phase transition curve (PTC) g(o(t.) has to be evaluated that tells neuron how its phase 6, changes when it receives an imput spike.
the phase representation, phases (t) evolve linearly in time between spikes, but the amount of change when the neuron receives an input spike depends on its phase. In this example at input spike time 11, the phase change is much larger than at input spike time 12, despite having same size in the voltage representation. In conventional event- based simulations, the phase all N neurons is shifted at every network spike time by Aé; WA,, where AL (ls+1 (.) is the time to the next network spike. Instead, in SparseProp just threshold and reset are changed at every spike.
through all neurons to find the next spiking neuron (step 1) and evolving all neurons to the next spike time (step 2) conventionally require O (N) calculations per network spike. In the following, we suggest an algorithm that for sparse networks only requires O (log(N) calculations per network spike by using change of reference frame and an efficient data structure. In the next paragraphs, we will describe these two features in more detail.
2.1 Change of Reference Frame
For clarity of presentation, we describe our algorithm first for a network of neurons that can be mapped pulse-coupled phase oscillators like the leaky integrate-and-fii neuron and the quadratic integrate- and-fire neuron (Fig [ A, C), and discuss later more general implementation. For pulse-coupled
oscillators, the network dynamics can be written as [10 15 16
17 18 19
f(x((t.) (t.) + w(s+11 - t.)) for i E post(j") (t.(t.) + w(s++) - )
(2)
else
where w is the constant phase velocity and g( is the phase transition curve (PTC) evaluated for the neurons that are postsynaptic to the spiking neuron j" just before the next network spike. In this phase representation, in the absence of recurrent input pulses, each phase evolves linearly in time with constant phase velocity from reset phase Greset to threshold phase 8th [ Thus, finding the next spiking neuron in the network j* amounts for homogeneous networks where all neurons have identical phase velocity to taking the maximum over all phases oz-
max(dx(t.))) (t.)).. The next Ph- g( quantifies the phase chånge neuron upon
network spike time 4++1 given by (y- function of its phase (Fig B). In this conventional form of the event- based algorithm [10 TO Ti 14] [5] 16, [7] [8] 19, 20 21). all N phases must be linearly evolved to the state immediately before the next network spike time t.+1; (t(++)) (t,)) w(t.-1 t,).
receiving an input spike
This operation, with computational cost of O (N), can be avoided in our proposed algorithm. Instead, we introduce a modification where just the threshold and reset phase are changed by a global phase offset Ag:
A (a+1) 1.)) + w(ttt1 - tt).
(3)
3

---

A
0.9
B
0.6
0.7
value
0.9
0.6
0.7
0.3
0.2
0.5
0.7
0.0
key
1
2
4
6
7
8
0.3
0.2
0.5
0.7
0.0
0.1
Figure 2: A suitable data structure allows efficient search for next spiking neuron for numeri- cally exact simulation of large spiking networks. By using priority queue, the next spiking neuron can be found without iteration through all neurons of the network at every spike time. One implementation of priority queue is binary heap. as the example shown. In a max-heap, each child node has value less or equal to its parent node. If an element is deleted, removed, or added, the heap property is restored by swapping parents and children systematically. Finding priority has numerical complexity of O (1) and changing the value of any element has amortized time complexity of (log(N) [)). Array implementation of binary heap shown in A. As per spike all postsynaptic neurons and the spiking neuron have to be updated, this requires on average K + operation, thus an amortized time complexity O (K log((N) 223
node highest
starting with A 0 at = 0. Thus 0 -
(s+1), Similarrly, Gaee greset (++1) See A
[i for illustration. In inhibition-dominated networks, the phases (t,,) will thus become increasingly negative. For long simulations, all phases and the global phase should therefore be reset if the global phase exceeds certain threshold Ath to avoid numerical errors resulting from subtractive cancellation due to Moating-poine arithmetic [22] We discuss this in more detail in appendix C In summary, this change of variables reduces the computational cost of propagating the network state from O (N) to O (1). Note, that assumed here identical phase velocity for all
neurons, but we will relax this assumption later section 6
2.2 Efficient Data Structure
We next address the computational bottleneck encountered in very large sparse networks when finding the next spiking neuron, which has a numerical time complexity of O (N). ). To overcome this limitation, we introduce an efficient implementation using priority queue. The priority queue is data structure that facilitates returning elements with the highest priority and supports insertion and priority updates 24 We implement the priority queue as heap-ordered binary tree, where each node's value is greater than its two children's 's values. This structure enables returning elements with the highest priority in constant time, with a time complexity of O (1), and inserting elements or updating priorities in logarithmic time, with time complexity of O (log(N)) [23] Thus, only O (K log( (N) operations have to be performed per network spike, as for postsynaptic neurons and for one spiking neuron, phase needs to be updated in the priority queue.

---

Algorithm SparseProp: Efficient Event-Based Simulation of Sparse Spiking Network
initialize o(to), 4.
2: heapify ((ta)
3: warm-up of network ((ta)
4: for s = = f
5: get phase of next spiking neuron: j. 6 (eck(
6: calculate phase increment: dø tt 0j (ta) + As
7:
update global phase shift: As += dø
8:
evaluate phase transition curve: 6 (t,) )= - ( (t,) + Ae)
9:
reset spiking neuron: ()(1++1) = 6e A
10:
if A > th then
11:
④ - A
12:
A. 0
13: end
14:
end for
where (te)( retrieves the highest-priority element without extracting from the heap. In utilized for the phase representation here -tee function fetches the smallest element, representing the next spike time. Conversely, a max-heap-employee it obtains the largest element, signifying the maximum phase. The operation generally executes in constant time, (1), given that the target element resides at the heap's root. Nonetheless, this function leaves the heap structure unaltered.
for heterogeneous networks-
An example implementation of SparseProp in Julia 25 is available here There are two core innovations in our suggested algorithm: First, we perform a change of variables into co-moving reference frame, which avoids iterating through all neurons in step 2 (Fig Second, for finding the next spiking neuron in the network, we put the network state in a priority heap (Fig A binary heap has efficient operations for returning the elements with the lowest (highest) key and insertion of new elements or updating of keys (24] We will next analyze the computational cost of SparseProp empirically and theoretically.
3 Computational Cost Scaling
Overall, the total amortized cost per network spike of SparseProp scales with O (K log(N)). For sparse networks, where K is fixed and only N is growing, the computational complexity of Sparse- Prop per network spike thus only scales with O (log(N)).
As the number of network spikes for given simulation time grows linearly with network size, the overall computational cost for given simulation time of SparseProp scales O (K N Iog(N)), which for sparse networks is far more efficient compared to a conventional implementation (see Fig. (3). Specifically, the linearithmic scaling with network size N of our approach, grants remarkable effi- ciency gains compared to the quadratic scaling of conventional event-based simulations, particularly for sparse networks. For example, simulating network of 10 neurons and K 100 synapses for 100s takes less then one CPU hour with SparseProp, but would take more then 3 CPU years with the conventional algorithn We provide detailed comparison of the computational cost in table
4 Example Implementation for Integrate-And-Firg Neuron Models
4.1 Pseudophase Representation of Leaky Integrate-And-Firg Neurons
For a network of pulse-coupled leaky integrate-and-In neurons, Eq. reads dimensionless notation
dV: = -V; ext +
Σ Jiy 8(t (()))
(4)
dl
3.5
'These numbers based on the benchmark Xeon@ CPU E5-4620 v2 @ 2.60 GHz and 512 GB RAM
5

---

conventional algorithm
SparseProp
Find next spiking neuron
min(óm @1/w)
peek(o:)
Evolve neurons
@ += dl
A+ Ae
Update postsynaptic neurons
K operations
K operations key updates O (log(N))
Reset spiking neuron
one array operation
O (log(N))
Memory cost
O (N)
O (N)
total amortized costs net- work spikk
O (N + K)
O (K log(())
total amortized costs for fixed simulation time
O (N²)
O (K N log(N))
Table 1: Comparison of computational cost for event-based simulation and training of recur- rent spiking networks for different algorithms N denotes number of neurons, K is the aver- age number of synapses per neuron. For large sparse networks (K < N) and a fixed number of synapses per neuron K, the dominant term grows quadratic for the conventional algorithm 28 and linearithmic (O log(N)) for the SparseProp algorithm.
10
26 27
B
108
100
cpu timm/siitlattin time
108
cpu time/spike (s)
104
10-2
102
conventional algorithm SparseProp algorithm
100
10-4
conventional algorithm SparseProp algorithm
10
10
10
10?
109
103
105
10
10
network size N
network size
Figure 3: Benchmark of SparseProp vs. conventional algorithm. The computational cost per network spike scales linear with network size N in a conventional event-based implementation. In the large sparse network limit, the computational bottleneck is to find the next spiking neuron in the network and to propagate all neurons phases to the next network spike. Both operations can be implemented efficiently. First, instead of shifting all neurons phases, the threshold and reset values can be shifted. Second, by using a binary heap as a data structure for the phases, finding the next phase and keeping the heap-ordering has computational complexity of O (K log(N)). A: CPU time per network spike as function of network size N for an inhibitory network of leaky integrate-and-ficg neurons with fixed number of synapses per neuron K [.. B: CPU time per simulation time shows quadratic scaling for conventional algorithm, but only linearithmic scaling for the novel Sparse agorithm, Benchmark was performed on an Intel@ Xeon@ CPU E5 E5-4
2.60 GHz and 512
GB RAM. Parameters: mean firing rate 7 Hz, Jo = 1, Tm 10,s, K = 100.
If a membrane potential V; reaches threshold Vnn it is reset to Via. Without loss of generality, we set V and V/2 Between two network spikes, the solution is given by:
vi(t+41) = rext (rext Vi((a)) exp (-
1.t1 (-1sta
fg
(5)
Tm
In this pseudophase representation (slightly different from [1 28 the phases di (-00, 0) describe the neuron states relative to the unperturbed interspike interval.

---

To obtain the unperturbed interspike interval 7 we have to solve Eq. in the absence of synaptic input.
between reset and threshold
Vín
Toxt
-Tm In
(6)
V/c
Text
Fm In
lext
Its inverse the phase velocity w = 1/T The phase di is thus given by Φ; = -w w
Jost
Text V;
The reverse transformation is V, = Text (1 - exp [- de ])
Therefore, the phase transition curve is
g (1++1)) = -w In (exp (-ox- (1--1)/w) +c)
where is the effective coupling strength c = Jext and the synaptic coupling strength taken here to be identical for connected neurons. Usually, as discussed later in the appendix Ejin the section on balanced networks, is scaled with J = Jn. example implementation of LIF network with SparseProp in Julia [25] is available here
VK' where K is the number of synapses per neuron. An
4.2 Phase Representation of Quadratic Integrate-And-Pi Neurons
The quadratic integrate-and- fire (QIF) neuron has, in contrast to the LIF neuron, dynamic spike generation mechanism and still has an analytical solution (i(t) between network spikes.
For network of QIF neurons, Eq. reads in dimensionless voltage representation:
dV Ti dt
V.2 + pext + Σ Jig 0(t J.8
(a)
(9)
The quadratic integrate-and-fir model can be mapped via a change of variables V (0/) to the theta model with a phase variable E (-T, (29 30 31. The dynamical equation between incoming spikes is the topological normal form for the saddle-node (SNIC) and allows closed-form solution of the next network spike thanks to the exact derivation of the phase response curve [32]. Therefore, the quadratic miegrate-and-fir neuron commonly used to analyze networks of spiking neurons [10] [3, 34 35 36). When 1 >> 0 A i, the right-hand
a limit cycle bifurcation
the absence of incoming postsynaptic potentials. In this case, we can choose another particularly tractable phase representation, called phi-representation with Vi(()= //ext tan (o: (())()))/)) where the neuron has constant phase velocity [10]. This transformation directly yields the phase transition curve
side of the dynamics is strictly positive and all neurons would spike periodically
2 Φi + +c).
g( arctan tan
(10)
where is the effective coupling strength d The phase velocity is given by Text, We considered homogeneous networks here, where all neurons have identical external input and identical membrane time constant, and will later consider the more general case of heterogeneous networks in section 6, An example implementation of QIF network with SparseProp in Julia (5 is available here
Event-Based Simulation via Chebyshev Polynomials
SparseProp can also be used for univariate neuron models where analytical solution between spike times is known, e.g., the exponential integrate-and-fin model (EIF) [77] 38 In this case, the phase transition curve can be calculated numerically before the network simulation. During run time, the phase transition curve and its derivative can be interpolated from the precomputed lookup tables or using Chebyshev polynomials (see appendix [B). An example implementation of an EIF network with SparseProp in Julia 25 is available here To test the numerical accuracy of the solution based on lookup table Chebyshev polynomials, we compared the spike times of two identical networks of LIF neurons with identical initial condition that were simulated either with the analytical
7

---

phase transition curve or using the numerically interpolated solution. We used LIF networks for this comparison, because pulse-coupled inhibitory LIF networks are non-chaotic [
39 40 4T]
Thus, numerical errors originating at the level of machine precision are not exponentially amplified unlike in the chaotic pulse-coupled QIF network dynamics [Q]. We show Fig. that the difference in spike time between the two networks remain close to machine precision and the temporal order
of the natwork enilae in altarad whan canlacino the onalutionl phoce troncition by the
1e-13
spike time enror Ar
2
spike index error
100000
100000
spike index
spike index
Figure 4: Left: Error of individual spike times for SparseProp are close machine precision. Right: Error in spike index. Note that despite the small enrors in the spike times that are close to machine precision, for non-chaotic network dynamics the spike index is still correct.
6 Time-based SparseProp for Heterogeneous Networks
For heterogeneous networks, where there is heterogeneity of external input current 1 or membrane time constant Tú, or threshold voltages Vig every neuron has a different phase velocity wi and single global phase shift described in section 2.1 does not suffice.
In that case, we suggest yet another representation. Instead of having a priority queue of phases, instead, we suggest priority queue of the next unperturbed spike time for all neurons, denoted by n(ta), this can take heterogeneity into account. this case, the algorithm for network simulation still consists of four steps: First finding the next spiking neuron, then updating the time followed by the update the postsynaptic neurons by the incoming spikes, finally, resetting the spiking neuron to the reset value. This algorithm also has numerical complexity of O (K log(N)) network spike. We provide pseudocode for SparseProp in the time-based representation below:
Algorithm 2 SparseProp: Event-Based Simulation of Heterogeneous Sparse Spiking Network
1: Initialize heterogeneous phase velocity w; based the neuron model
2: Initialize unperturbed interspike interval T;
3: Initialize time to next spike for all neurons n(to)
4: Initialize global time shift A, =
5: Heapify n(fo)
6: Perform warm-up of network state n(to)
7: for 1
8: Get index and spike time of next spiking neuron: j. ny peek(n((.))
9: Calculate time increment: dt (1y(1s) A,
10: Update global time shift: A, += dt
Update postsynaptic neurons i": n,t (t.) update (n, (to), A.)
12:
Reset spiking neuron: (y(t++1) A TV
13:
if A 4 then
14:
n A:
15:
A,
16:
end
17:
end for
Again, (ee)(n((.)) retrieves the highest-priority element without extracting it from the heap, which in this case is the minimum over all neurons' next unperturbed spike time. An example implementation of heterogeneous QIF network with SparseProp in this time-based representation is available here

---

Note that this can also easily be extended to mixed networks of excitatory and inhibitory neurons (or population networks with different external input 42
The update of the postsynaptic neurons in line 11 is analogous to the evaluation of the phase transition curve in the homogeneous SparseProp, but takes into account the different phase velocities. For example, in the case of QIF neuron, it takes the form
π g(TT (n: - A A))we) + A.
apdateenn, (t. (t.).
(11)
Wi
Here, the phase transition curve, denoted as g( is the same as previously described in Eq. 10 g(or) arctan (tan &e + c) with effective coupling strength For the QIF neuron, the period T; is given by T, 2T/ws, with w = 2
7 Efficient Event-Based Training of Spiking Neural Networks
SparseProp can also be used for efficient event- based training of spiking neural networks. While
46), it was recently claimed that EventProp, which uses exact analytical gradients based on spike times, can also yield competitive performance [47]. To obtain precise gradients for event-based simulations, this algorithm makes use of the adjoint method from optimization theory 48 49]. In this approach, in the forward pass the membrane potential dynamics (Eq. I is solved in event-based simulation and spike times are stored. In the backward pass, an adjoint system variables is integrated backward in time, and errors in spike times back propaprop through the size of the jumps in these adjoint variables simulations of neuron models where the next spike time
surrogate gradients seem to facilitate training spiking networks 43 44
45
at the stored spike times. For is analytically known or be obtained like in section [5] the entries adjoint integration is being taken care of by differential programming [47]. Besides the previously mentioned advantages of event-based simulation regarding precision, this approach is not only sparse in time but also sparse with respect to the synapses, as only synaptic events are used.
We suggest here further improve the efficiency of event-based training using SparseProp. As the backward pass has the identical amortized time complexity of O [log(N), per network spike as the forward pass, we expect significant boost of training speed for large sparse networks 447 [0 [T]. Regrettably, the implementation of EventProp presented in [4]] was not publicly available during the time of submission. In contrast, the authors of (50) made their re-implementation of EventProp publicly available, but they did not use an event- event-based implementation but forward Euler to integrate the network dynamics. It will be important next step to also benchmark the performance of
SparseProp training. As the gradients are identical, we expect similar results to 47
50 51
Limitations
In this work, focus exclusively recurrent networks of pulse-coupled univariate neuron models such as the quadratic integrate-and-fin neuron, Ieaky integrate-and fire neuron, and the exponential integrate-and and-fire neuron thereby excluding the exploration of multivariate neuron models. Extending our approach to encompass multivariate models, such as the leaky- or quadratic integrate-and -ire remains a promising direction for future research.
neurons with slow synapses [99 28
51
more fundamental limitation of our study is the apparent incompatibility of event- based simulations with surrogate gradient techniques 466 53 While could introduce 'ghost spikes' in the event- based simulation to emulate surrogate gradients when neurons hit a lower threshold. However, it remains unclear how preserve the favorable computational scaling of SparseProp in this case.
Furthermore, SparseProp requires an analytical solution of Vt(t) given the initial state and the neuronal dynamics (1()) is necessary. This limitation excludes conductance- based models [54] 55 (56. for example. Additionally, time- varying input only be introduced in the form of deterministic currents or point processes, such as Poisson input. Incorporating stochastic input without an analytical solution of Vt(l) would require further development.
Moreover, the reduction in computational cost per network spike from (N) to O (log(N) can only be achieved in sparse networks, where the number of synapses per neuron K is much smaller than the

---

total number of neurons N. For dense networks where the number of synapses scales proportionally
to the number of neurons [B] a priority queue implemented by a binary heap is disadvantageous compared to conventional array in the large network limit, as every network spike involves changing the priority of O (N) neurons thus O ((N + 1) log(N) flops per network spike which corresponds
to O ((N + 1) N log (og( flops for fixed simulation time. A batch-update of all postsynaptic neurons might be faster for very large networks [57] but this is beyond the scope of this work. would involve O (K + log(K) K) log(N)) flops corresponding to O (N + log(N) log(N)) in dense networks. In the case of purely excitatory sparse spiking networks, a Fibonacci heap might be a more efficient implementation in the large network limit [58], the 'decrease key' operation takes constant time O (1) compared to O (l) in the case of a binary heap. Note that for practical purposes, the asymptotic scaling of the computational complexity of Fibonacci heaps has unfavorably large prefactor [59]. Therefore, other heap structure implementations might be faster [60]. While large mammalian brains sparse N > K, for models of small brains or local circuits with dense connectivity the conventional event-based simulations might be faster.
Additionally, SparseProp cannot currently model biophysically detailed neuron features, including dendrites, threshold adaptation, slow synaptic timescales, and short-term plasticity. Incorporating these features would be another valuable avenue for future computational neuroscience research.
Lastly, we did not yet perform a comprehensive systematic benchmark of training with challenging tasks. Such benchmarks would offer further insights into the training performance of our proposed framework.
Discussion
We introduce SparseProp, an efficient and numerically exact algorithm for simulating and training
large sparse spiking neural networks in event-based simulations. By exploiting network sparsity, we achieve a significant reduction in computational cost per network spike from O (N) to O [log(N)) This speedup achieved by optimizing two critical steps in the event-based simulation: finding the next spiking neuron and evolving the network state to the next network spike time. First, we employ a binary heap data structure for finding the next spiking neuron and updating the membrane potential of postsynaptic neurons. Second, we utilize a change of variables into co-moving reference frame, avoiding the need to iterate through all neurons at every spike.
Our results demonstrate the scalability and utility of SparseProp in numerical experiments. We demonstrate that SparseProp speeds up the simulation and training of a sparse SNN with one million neurons by over four orders of magnitude compared to previous implementations. In contrast to the conventional event-based algorithm, which would require three CPU years to simulate 100 seconds of network dynamics, SparseProp achieves the same simulation in just one CPU hour.
This advancement enables the simulation and training of more biologically plausible. large-scale neural network models relevant to theoretical neuroscience. Furthermore, it might pave the way for the exploration of segle-scaleeae event- spiking neural networks in other machine learning domains, such as natural language processing [61]. embedded automotive applications
622 robotics (63] and
in silico pretraining of neuromorphic hardware [6.4.
64].
The impact of SparseProp extends to the recent surge of training for spiking neural networks in theoretical neuroscience [50) 65), machine learning [61, 666 67). and neuromorphic computing
668 69 70 This algorithm is expected to substantially accelcrate the training speed in recent works focused on event-based spiking networks
50 51].
Future extensions of SparseProp to multilayer and feedforward networks would be a promising avenue. Integration of SparseProp with recent attempts in biologically more plausible backpropa- gation and real-time recurrent learning [53] 65
72 74 could be used in future research in theoretical neuroscience. To improve training performance, future studies should consider incorporat- ing additional slow time scales, such as slowly decaying synaptic currents 722, and slow threshold adaptation (65), which appear to enhance temporal-credit assignment in recurrent neural networks for rudimentary tasks [53 65 72). Additionally, incorporating gradient shaping and insights from advances in event-based spiking network training [50] [1) could yield further gains in performance.
10

---

Acknowledgments and Disclosure of Funding
I thank S. Goedeke, J. Liedtke, R.-M. Memmesheimer, M. Monteforte, A. Palmigiano, M. Puelma Touzel, A. Schmidt, M. Schottdorf, and F. Wolf for fruitful discussions. Moreover, I thank L.F. Abbott, Migirditch,
Palmigiano, Y. Park, and D. Sussillo for comments on the manuscript. Finally, thank the Julia community.
Research supported by NSF NeuroNex Award (DBI-1707398), the Gatsby Charitable Foundation (GAT3708)) the Simons Collaboration for the Global Brain (542939SPI), and the Swartz Foundation (2021-6). This work further supported by the Deutsche Forschungsgemeinschaff (DFG, German Research Foundation) 436260547 in relation to NeuroNex (National Science Foundation 2015276) and under Germany's Excellence Strategy EXC 2067/1- 390729940, DFG Project-ID 317475864 SFB 1286, by Evangelisches Studienwerk Villigst, by the Leibniz Association (project K265/2019), and by the Niedersächsisches Vorab of the VolkswagenStiftung through the Göttingen Campus Institute for Dynamics of Biological Networks.
References
[1] Peter U. Diehl, Guido Zarrella, Andrew Cassidy, Bruno U. Pedroni, and Emre Neftci. Conversion of artificial recurrent neural networks to spiking neural networks for low-power neuromorphic hardware. 2016 IEEE International Conference on Rebooting Computing (ICRC), pages 1-8, October 2016.
[2] Payal Dhar. The carbon impact of artificial intelligence. Nature Machine Intelligence, 2(8):423- 425, August 2020.
[3] Tim Gollisch and Markus Meister. Rapid Neural Coding in the Retina with Relative Spike Latencies. Science, 319(5666): 108-1111,,,11 February 2008.
[4] Robert Gütig and Haim Sompolinsky. The tempotron: neuron that learns spike timing-based based decisions. Nature Neuroscience, 9(3):420-444 428, March 2006.
[5] Robert Gütig. Spiking neurons can discover predictive features by aggregate-label learning. Science, 351(6277):ab44133 March 2016.
[6] D. Hansel, G. Mato, C. Meunier, and L. Neltner. On Numerical Simulations of Integrate-and-Fire Neural Networks. Neural Computation, 10(2):467-4833 February 1998.
[7] D. J. Amit, M. R. Evans, and M. Abeles. Autractor neural networks with biological probe records. Network: Computation in Neural Systems, 1(4):811, October 1990.
[8] M V Tsodyks and T Sejnowski. Rapid state switching in balanced cortical network models. Network: Computation in Neural Systems, 6(2):1 January 1995.
[9] Nicolas Brunel. Dynamics of Sparsely Connected Networks of Excitatory and Inhibitory Spiking Neurons. Journal of Computational Neuroscience, 803):183-208, May 2000.
[10] Michael Monteforte and Fred Wolf. Dynamical Entropy Production in Spiking Neuron Networks in the Balanced State. Physical Review Letters, 105(26):268 104, December 2010.
[11] Michael Monteforte and Fred Wolf. Dynamic Flux Tubes Form Reservoirs of Stability in Neuronal Circuits. Physical Review X, 2(4):041007, November 2012.
[12] Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski. Neuronal Dynamies: From Single Neurons to Networks and Models of Cognition. Cambridge University Press, Cambridge, uk ed. edition edition, September 2014.
[13] Alfonso Renart, Jaime de la Rocha, Peter Bartho, Liad Hollender, Néstor Parga, Alex Reyes, and Kenneth D. Harris. The Asynchronous State in Cortical Circuits. Science, :9755965),8775590 January 2010.
[14] Agostina Palmigiano, Rainer Engelken, and Fred Wolf. Boosting of neural circuit chaos the onset of collective oscillations. eLife, 12, November 2023.
11

---

[15] Renato E. Mirollo and Steven H. Strogatz. Synchronization of Pulse-Coupled Biological Oscillators. SIAM Journal on Applied Mathematics, 5066):6645-162,, December 1990.
[16] M. Tsodyks, I. Mitkov, and H. Sompolinsky. Pattern of synchrony in inhomogeneous networks of oscillators with pulse interactions. Physical Review Letters, 1788))2200-1288 August 1993.
[17] U. Ernst, Pawelzik, and T. Geisel. Synchronization Induced by Temporal Delays in Pulse- Coupled Oscillators. Physical Review Letters, 74(9):1570-1573, February 1995.
[18] Marc Timme, Fred Wolf, and Theo Geisel. Unstable attractors induce perpetual synchronization and desynchronization. Chaos (Woodbury, N.Y.), 13(1::37-387, March 2003.
[19] Michael Monteforte. Chaotic Dynamics in Networks of Spiking Neurons in the Balanced State. PhD thesis, Boong-August-Univeee Göttingen, May 2011.
[20] Alexander Schmidt. Characterization of the phase space structure in driven neural networks in Göttingen, 2016.
the balanced state. Master's thesis, Lrrrg--Aguu-t-nnn
[21] Maximilian Puelma Touzel and Fred Wolf. Statistical mechanics of spike events underlying phase space partitioning and sequence codes in large-scale models of neural circuits. Physical Review E, 9955)0052 May 2019.
[22] David Goldberg. What every computer scientist should know about floating-por -pinnt arithmetic. ACM Computing Surveys (CSUR). 23(1):5-48, March 1991.
[23] Robert Sedgewick and Kevin Wayne. Algorithms. Addison- Addison-Wesley Professional, Upper Saddle River, NJ, 4th edition edition, March 2011.
[24] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. Introduction to Algorithms, 3rd Edition. The MIT Press, Cambridge, Mass, 3rd edition edition, July 2009.
[25] J. Bezanson, A. Edelman, S. Karpinski, and V. Shah. Julia: A Fresh Approach to Numerical Computing. SIAM Review, 59(1):65-99) 98, January 2017.
[26] Dylan Festa. Chaos Characterization of Pulse-Coupled Neural Networks in Balanced State. Master's thesis, MPI DS / Università di Pisa, Göttingen/Pisa, 2011.
[27] Joscha Liedtke. Geometry and organization of stable and unstable manifold in balanced networks. Master's thesis, Georg-August- Augusss-U----vsessttt Göttingen, 2013.
[28] Maximilian Puelma Touzel. Cellular dynamics and stable chaos in balanced networks. PhD thesis, Aooog--uuust-Universtt Göttingen, January 2016.
[29] G. Ermentrout and N. Kopell. Parabolic Bursting in an Excitable System Coupled with Slow Oscillation. SIAM Journal on Applied Mathematics, 46(2):233-253, April 1986.
[30] Bard Ermentrout. Type I Membranes, Phase Resetting Curves, and Synchrony. Neural Computation, 8(5):979-1001, July 1996.
[31] B. S. Gutkin and G. B. Ermentrout. Dynamics of Membrane Excitability Determine Interspike Interval Variability: A Link Between Spike Generation Mechanisms and Cortical Spike Train Statistics. Neural Computation, 1065):1047-1055, July 1998.
[32] Eugene M. Izhikevich. Dynamical Systems in Neuroscience. MIT Press, 2007.
[33] Remus Osan and Bard Ermentrout. Two dimensional synaptically generated traveling waves in
a theta-neuron neural network. Neurocomputing, 38 40:789-79$ June 2001.
[34] B. S. Gutkin, J. Jost, and H. C. Tuckwell. Transient termination of spiking by noise in coupled neurons. EPL (Europhysics Letters), (1)2:200055 2008
[35] D. Hansel and G. Mato. Asynchronous States and the Emergence of Synchrony in Large Networks of Interacting Excitatory and Inhibitory Neurons. Neural Computation, 15(1):1-56. January 2003.
12

---

[36] P. E. Latham, B. J. Richmond, S. Nirenberg, and P. G. Nelson. Intrinsic Dynamics in Neuronal Networks. II. Experiment. Journal of Neurophysiology, 83(2):828-855 February 2000.
[37] Nicolas Foureaud- Troemé, David Hansel, Carl van Vreeswijk, and Nicolas Brunel. How Spike Generation Mechanisms Determine the Neuronal Response to Fluctuating Inputs. The Journal of Neuroscience, 23(37): 1668)))11228--1 11640, December 2003.
[38] Romain Brette. Exact Simulation of Integrate-and Models with Exponential Currents. Neural Computation, 19(10):2004-260) August 2007.
[39] Sven Jahnke, Raoul-Martin Memmesheimer, and Mare Timme. Stable Irregular Dynamics in Complex Neural Networks. Physical Review Letters, 100(4):048102, January 2008
[40] Sven Jahnke, Raoul-Martin Memmesheimer, Mare Timme, Sven Jahnke, Raoul- Memmesheimer, and Mare Timme. How chaotic is the balanced state? Frontiers in Computa- tional Neuroxcience, 3:13, 2009.
[41] Rüdiger Zillmer, Nicolas Brunel, and David Hansel. Very long transients, irregular firing, and chaotic dynamics in networks of randomly connected inhibitory integrate-and- and-fire neurons. Physical Review E, 79(3):031909, March 2009.
[42] Tobias C. Potjans and Markus Diesmann. The Cype-Type Specific Cortical Microcircuit: Relating Structure and Activity in Full-Scale Spiking Network Model. Cerebral Cortex, December 2012.
[43] Sander M. Bohte. Error- Bane-Baacpropaatioon Networks of Fractionally Predictive Spiking Neurons. In Timo Honkela, Włodzislaw Duch, Mark Girolami, and Samuel Kaski, editors, Artificial Neural Networks and Machine Learning ICANN 2011, Lecture Notes in Computer Science, pages 60-68, Berlin, Heidelberg, 2011. Springer.
[44] Sander M. Bohte, Joost N. Kok, and Han La Poutré. Error-backproppggtion in temporally encoded networks of spiking neurons. Neurocomputing. 48(1):17-37, October 2002.
[45] Olaf Booij and Hieu tat Nguyen. A gradient descent rule for spiking neurons emitting multiple spikes. Information Processing Letters, 95(6):552-558, September 2005.
[46] Emre O. Neftci, Hesham Mostafa, and Friedemann Zenke. Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks. IEEE Signal Processing Magazine, 36(6):51-633 2019.
[47] Timo C. Wunderlich and Christian Pehle. Event-based backpropagation can compute exact gradients for spiking neural networks. Scientific Reports, 11(1):12829, June 2021.
[48] L. S. Pontryagin. Mathematical Theory of Optimal Processes: The Mathematical Theory of Optimal Processes. Routledge, New York, 1st edition edition, March 1987.
[49] Daniel Liberzon. Calculus of Variations and Optimal Control Theory: A Concise Introduction. In Calculus of Variations and Optimal Control Theory. Princeton University Press, December 2011.
[50] Thomas Nowotny, James P. Turner, and James C. Knight. Loss shaping enhances exact gradient learning with EventProp in Spiking Neural Networks. Technical report, December 2022. ADS Bibcode: 2022arXiv22 120 Type: article.
[51] Christian Klos and Raoul-Martin Memmesheimer. Smooth Exact Gradient Descent Learning in Spiking Neural Networks. 2023.
[52] Arnaud Tonnelier, Hana Belmabrouk, and Dominique Martinez. Event-Driven Simulations of Nonlinear Integrate-and-Fing Neurons. Neural Computation, 19(12):326-33238 October 2007.
[53] Friedemann Zenke and Surya Ganguli. SuperSpike: Supervised Learning in Multilayer Spiking Neural Networks. Neural Computation, 30(6): June 2018
[54] A. L.. Hodgkin and A. F. Huxley. The components of membrane conductance in the giant axon of Loligo. The Journal of Physiology, 1 116(4):433-96 April 1952.
13

---

[55] A. L. Hodgkin and A. F. Huxley. Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo. The Journal of Physiology, 116(4):449-47)) April 1952.
[56] Xiao-Jing Wang and György Buzsáki. Gamma Oscillation by Synaptic Inhibition in a Hip- pocampal Interneuronal Network Model. The Journal of Neuroscience, 16(20):6402-6413 October 1996.
[57] Jörg-R. R. Sack and Thomas Strothotte. An algorithm for merging meaps. Acta Informatica. 22(2):171-186, June 1985.
[58] Michael L. Fredman and Robert Endre Tarjan. Fibonacei heaps and their uses improved network optimization algorithms. Journal of the ACM (JACM), 34(3):596-615, July 1987.
[59] Michael L. Fredman, Robert Sedgewick, Daniel D. Sleator, and Robert E. Tarjan. The pairing heap: new form of self-adjusting heap. Algorithmica, 1(1-4):111-13 November 1986.
[60] Gerth Stølting Brodal. Worst-case efficient priority queues. Proceedings of the Seventh Annual Acm-siam Symposium on Discrete Algorithms, pages 52-58, 1996.
[61] Rui-Jie Zhu, Qihang Zhao, and Jason K. Eshraghian. SpikeGPT: Generative Pre-trained Language Model with Spiking Neural Networks. Technical report, February 2023. ADS Bibcode: 2023arXiv2302 3939Z Type: article.
[62] Loie Cordone, Benoît Miramond, and Philippe Thierion. Object Detection with Spiking Neural Networks on Automotive Event Data. In 2022 International Joint Conference on Neural Networks (IJCNN) pages 1-8, July 2022. ISSN: 2161-4407.
[63] Travis Wolf Spiking neural networks take control. Science Robotics, 6(58))eabk3368. September 2021.
[64] Eric Hunsberger and Chris Eliasmith. Training Spiking Deep Networks for Neuromorphic Hardware. 2016. arXiv:161 1.056611111144 141 [cs].
[65] Guillaume Bellec, Franz Scherr, Anand Subramoney, Elias Hajek, Darjan Salaj, Robert Leg- enstein, and Wolfgang Maass. A solution to the learning dilemma for recurrent networks of spiking neurons. Nature Communications, 11(1):3625, July 2020.
[66] Iulia-Maria Comşa, Krzysztof Potempa, Luca Versari, Thomas Fischbacher, Andrea Gesmundo, and Jyrki Alakuijala. Temporal Coding in Spiking Neural Networks With Alpha Synaptic Function:
Learning With Backpropagation. IEFE Transactions Neural Networks and Learning Systems, (3(10))5939-599 5552,, October 2022.
[67] Bojian Yin, Federico Corradi, and Sander M. Bohté, Accurate online training of dynamical spik- ing neural networks through Forward Propagation Through Time. Nature Machine Intelligence, pages 1-10. May 2023.
[68] Paul A. Merolla, John V. Arthur, Rodrigo Alvarez-lcaza, Andrew S. Cassidy, Jun Sawada, Filipp Akopyan, Bryan L. Jackson, Nabil Imam, Chen Guo, Yutaka Nakamura, Bernard Brezzo, Ivan Vo, Steven K. Esser, Rathinakumar Appuswamy, Brian Taba, Arnon Amir, Myron D. Flickner, William P. Risk, Rajit Manohar, and Dharmendra S. Modha. A million spiking-neuron integrated circuit with a scalable communication network and interface. Science, 345(6197):668- 673, August 2014.
[69] J. Göltz, L. Kriener, A. Baumbach, S. Billaudelle, O. Breitwieser, B. Cramer, D. Dold, A. F. Kungl, W. Senn, J. Schemmel, K. Meier, and M. A. Petrovici. Fast and energy-efficient neuromorphic deep learning with first-spike times. Nature Machine Intelligence, 32))8823-83, September 2021.
[70] Christian Pehle, Luca Blessing, Elias Arnold, Eric Müller, and Johannes Schemmel. Event- based Backpropagation for Analog Neuromorphic Hardware. Technical report, February 2023. Type: article.
ADS Bibcode: 2023ariiw2302074411
[71] Owen Marschall, Kyunghyun Cho, and Cristina Savin. A unified framework of online learning algorithms for training recurrent neural networks. The Journal of Machine Learning Research, 2I(D):135:3200-135555 135:5333, January 2020.
14

---

[72] James M Murray, Local online learning in recurrent networks with random feedback. eLife, 8:e43299, May 2019.
[73] Alexandre Payeur, Jordan Guerguiev, Friedemann Zenke, Blake A. Richards, and Richard Naud. Burst-dependent synaptic plasticity can coordinate learning in hierarchical circuits. Nature Neuroscience, 24(7):1010-1019) July 2021.
[74] Yuhan Helena Liu, Stephen Smith, Stefan Mihalas, Eric Todd SheaBrown, and Uygar Sümbül. moologically-plausiie backpropagation through arbitrary timespans via local neuromodulators. October 2022.
[75] Paul Manz, Sven Goedeke, and Raoul-Martin Memmesheimer. Dynamics and computation in mixed networks containing neurons that accelerate towards spiking. Physical Review. E, 100(4-1);042404, October 2019.
[76] Eugene M. Izhikevich and Gerald M. Edelman. Large-scale model of mammalian thalamocorti- cal systems. Proceedings of the National Academy of Sciences, 105(9):3593-3598 2008
[77] Robert Rosenbaum and Brent Doiron. Balanced Networks Spiking Neurons with Spatially Dependent Recurrent Connections. Physical Review X, 4(2):021039, May 2014.
[78] William R. Softky and Christof Koch. Cortical Cells Should Fire Regularly, But Do Not. Neural Computation, 4(5):643 646, September 1992.
[79] W. R. Softky and C. Koch. The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. The Journal of Neuroscience, 13(1):334 350, January 1993.
[80] W. H. Calvin and C. Stevens. Synaptic noise and other sources of randomness in motoneuron interspike intervals. Journal of Neurophysiology, 31(4):574-587, July 1968.
[81] H Bryant and P Segundo. Spike initiation by transmembrane current: a white-noise analysis. The Journal of Physiology, 260(2):27933) September 1976.
[82] John D. Hunter, John G. Milton, Peter J. Thomas, and Jack D. Cowan. Resonance Effect for Neural Spike Time Reliability. Journal of Neurophysiology. 8033):1277-1438, September 1998.
[83] Z. F. Mainen and T. J. Sejnowski. Reliability of spike timing in neocortical neurons. Science, (5055226):50031550 June 1995.
[84] Michael N. Shadlen and William T. Newsome, Noise, neural codes and cortical organization. Current Opinion in Neurobiology, 4(4):569-579, August 1994.
[85] Michael N. Shadlen and William T. Newsome. The Variable Discharge of Cortical Neurons: Im- plications for Connectivity, Computation, and Information Coding. The Journal of Neuroscience, 18(10):380-3386 May 1998.
[86] C. van Vreeswijk and H. Sompolinsky. Chaos in Neuronal Networks with Balanced Excitatory and Inhibitory Activity. Science, 274(5293):1724 -1726, December 1996.
[87] C. van Vreeswijk and H. Sompolinsky. Chaotic Balanced State in a Model of Cortical Circuits. Neural Computation, 1066):1321-1371 1998.
[88] Tom Tetzlaff, Moritz Helias, Gaute T. Einevoll, and Markus Diesmann. Decorrelation of Activity by Inhibitory Feedback. PLOS Comput Biol, 1002599, August
2012.
15

---

Appendix
A SparseProp Code for Efficient Spiking Network Simulation in Julia
This minimal code implements the SparseProp algorithm for simulating spiking network of N = -00 LIF neurons.
using DataStructures, RandomNumbers: Xorshifts, StatsBase, PyPlot
function lifnet (n, nsteepp,, jo, r,sedii, ,eeddic,
iext r*sqrt (k) j00
iext given by balance equation phase velocity LIF
w,c (1 1/iext), j0/sqrt (k)/ (1. + iext)
øth, óshift 1., 0.
threshold for LIF
I Xoroshiro128Plus (seedic)
# init. random number generator
# initialize binary heap:
á Matabll@inaryJieaayPPooott
DataStructures. Res.:aurreerrtttteteeressssse (rand(r,n))
spikeidx = Int64 []
#initialize time # spike raster
spiketimes Float64
postidx rand (Int,k
for nstep
# main loop
omax, top_with_handleede
get phase next spiking neuron # calculate next spike time global shift evolve network state # spiking neuron index is seed of rng # get receiving neuron index
dø = Øth çmax @shift
øshift
Random. seed! ! (r, (r,jesee)net) sample! (r, 1:n-1, ptt.iod p
@inbounds for i postidx >= & ( [ostidx[i]++s]] )
# avoid autapses
end ptc! poostixxddshitt,..o (O.
# evaluate phase transition curve reset spiking neuron
push! (spiketimes, øshift) push! (spikeidz,j)
# store spike times
# store spiking neuron index
end
# output: rate, spike times indices
end
function ptc! (Ø, postid, øshift, w, c)
# phase transition curve of LIF
for postid
end [[i] w+log(exp( (d[i) dshift)/w) + c) øshift #(Eq. 12)
end
set parameters:
#n: # of neurons, k: synapses/neuron, j0: syn. strength, т: membr. time const. n, jo,, trratesnt,,,seddic,,,t seednet 10"5,10'0 5, 100, 1, 11, quick run to compile code
@time lifnet 1, 10, jo, ratewnt, т, seedic, seednet);
run benchmark network with specified parameters GC. (e;.gc();timm rate, sidx, stimes = lifnet (n, nsteptep k,
plot spike raster plotolstiii (stimes sitimms,,idx,, ,,k" ms=0..) ylabel ("Neuron Index" ,fontsize=20)
xlabel ("Time fontssiieeee=))) layout ()
16

---

B Event-Based SparseProp With Numerical Phase Transition Curve
For neuron models where the next spike time cannot be determined analytically and/or the net- work state cannot be temporally evolved based on a closed-form expression V(ta+z), such as the exponential integrate-and-
model, we present an event-based implementation that has the same computational cost of O (K log(N)) per network spike. To achieve this, we evaluate the phase transition curve either by employing a lookup table or by approximating it using Chebyshev poly- nomials. We initially describe the implementation based on a lookup table and provide an example implementation:
We begin by numerically integrating the ordinary differential equation that describes the membrane potential V of the neuron model, without considering any external input spikes:
dV
Tm dt
F(V)
text.
(12)
The integration is performed from the reset voltage Vie to the threshold voltage Vn as follows:
V(I) = I. (P(V) + )/m-
(13)
We obtain the numerical solution (t)) by employing the DifferentiallEuatiin.. j1 package with a callback that terminates the integration precisely at Va-
The solution provides us with two crucial pieces of information. Firstly, for each time point I, we can now associate a corresponding voltage (I). Secondly, we also acquire the unperturbed interspike interval T, which determines the phase velocity w. Using the numerically obtained solution ((), we create a lookup table that maps time to voltage by employing the Data.nterpolation package. By taking the inverse, we obtain mapping from a voltage to the next spike time t(V), which we also store as a separate lookup table. The two lookup tables combined give us a numerical phase transition curve g We contirmed for the leaky integrate-and-fin neuron and the quadratic integrate-and and-ffre neuron that this numerical procedure indeed yields the correct phase response curve d. Subsequently, incorporate the numerical phase response curve d into the event-based implementation SparseProp implementation and also confirmed the spike times remained identical when comparing with LIF networks with the analytical phase transition curve and identical network realization and identical initial conditions. We used an inhibitory LIF network for this sanity check, because inhibitory pulse-coupled networks of leaky integrate-and-fin neuron stable with respect to infinitesimal perturbations [M] 22 39 40 41]. This feature in combination with instability with respect finite-size perturbations renders inhibitory LIF networks an ideal test-based for numerical precision, because event-based simulations should result exactly in the same spike sequence. This is the case for networks of QIF neurons, which typically chaotic [10 (see also [75
Algorithm 3 SparseProp:: Eivent-Based Simulation Based on Lookup Tables
1: Set up the ordinary differential equation Eq. 12
2: Numerically integrate the ODE from Vr= to Vin to obtain V(I).
3: Generate a lookup table based on the obtained solution V()).
4: Generate inverse lookup table based on the solution t(V).
5: Combine the lookup table and inverse lookup table to obtain the phase transition curve g()).
6: Heapify o(to)
7: Perform the warm-up of network ((o))
8: for + do
9: Get index and phase of next spiking neuron: d. = (te)(
10:
Calculate phase increment: dó 6t + - (t.)
11:
Update global phase shift: A += dg
12:
Evaluate phase transition curve: o (t,) = ( (t.)) + A)
13:
Reset spiking neuron: (t(+++1) Th A
14:
if
> A th then
15:
+ A A
16:
17:
end if
18:
end for
17

---

An example implementation of an EIF network using SparseProp in Julia, in the balanced state, is available here
As an alternative to generating two lookup tables to numerically define the phase transition curve g using DataInterpolationa j1,, propose the use of Chebyshev Polynomials for direct estimation of d. We provide an example implementation using ApproxFun. here
C Numerical precision of SparseProp
Since SparseProp employs event-based simulation with analytical expressions for calculating the next spike times, the algorithm operates with a precision limited only by machine epsilon. Numerically, over many spike events, the global phase offset Aø incrementally increases. For sufficiently large values of Δć, catastrophic cancellation can arise. For instance, when Áó 10t and this is subtracted from very small dø. Under these conditions, the issue is that the adjacent floating-poit greater than 10t in double-precistee
number
arithmetic differs by approximately 2 -30 ze 1.16 x 10 10, One can mitigate this numerical error by resetting both the global phase shift and the phases individual neurons upon exceeding specified threshold. This counteracts the subtractive cancellation errors inherent floating- -ointig calculations. Nonetheless, such resets are seldom required and occur at intervals on the order of cN steps. The coefficient e is substantial prefactor contingent on the tolerable floating-poin error; for an error of -30, the command log2 (nextfloat (1. 0e6) -1 0e6).
-30, e would 10 This can be verified Julia with
D Online Sparse Erdős-Rényi Connectivity Matrix Generation
order to benchmark the simulation of large networks using the SparseProp algorithm Fig we implemented online generation of the sparse Erdös- connectivity matrix. This approach allows us to generate the set of postsynaptic neurons deterministically online, without the need to store the entire weight matrix. Storing the entire adjacency matrix for networks with 10 neurons would be impractical due to its prohibitively large size, even when using the compressed sparse column format. By using this online generation method, only require memory to store the state of all N neurons. For a network 10 neurons, this corresponds to approximately 64bit/double float x ne x 8bytes/bit 7.4 GB of RAM. Our implementation draws inspiration from previous works [76] (77). is important to note that the example code provided in appendix A has a fixed indegree, similar to the approach in 9 Here, present example of online sparse adjacency matrix generation with truly Frdös-Rényí connectivity, where the indegree follows binomial distribution.
Code example:
18

---

This code for Julia v1.9 implements sparse Erdós-Rény: Rényi connectivity matrix for efficiently simulat- ing a spiking network.
using Random, DataStructures, RandomWumbers.. Xorshifts, StatsBase
# before simulation:
N, K, J, seedlet 10"6, 10"2, 2,
init. random number generator (Xorshift) npostAll Array(UInt11,, (undef, N)
rng (rr
for n 1:N upostAll [n] UInt16 (raeddirrnnbin N, K N))
end
before simulation when neuron j spikes: npost upostAll[j]
spiking neuron index is seed of rng Random seed! (rng, j)
get postsynaptic neuron index postIdxes sample(rng, 1:N-1, npost, replace false)
# avoid autapses
for np 1:npost
if postIdxes [np] >=
postIdxes [np]
end
end
before simulation, the following function was used function randbinom (rng, n: n: Integer, p::Real)
log_q = log(1.0 p)
sum 0.0
while true
sum += log (rng) (n x)
sum log_q & break x +=
end
return x
end
E
Spiking Networks in the Balanced State
In order to provide a comprehensive understanding and ensure the accessibility of our work, we now introduce the concept of the balanced state, which motivates the scaling of the synaptic coupling strength with 1/VK. Neural activity cortical tissue has typically asynchronous and irregular pattern of action potentials [78 77). despite the fact that individual neurons can respond reliably This is commonly explained by a balance of excitatory and inhibitory synaptic which cancels large mean synaptic inputs average. A dynamically self-organized balance can be achieved without the fine-tuning of synaptic coupling strength in heterogencous networks, if the connectivity is ciaibition-domnnatee and the couplings are strong such that small active fraction of incoming synapses can trigger an action potential [86]. The statistical properties of this state are described by mean-field theory, which is largely insensitive to the specific neuron model (87.
80 currents
81 82 83
In our study, we investigated large sparse networks consisting of N neurons arranged on directed Frdōs-Rény: random graph with mean degree of K. All neurons 1,..., received constant positive externa currents eext and non-delayed d-pulses from the presynaptic neurons pre (i).
19

---

The values of the external currents were chosen using bisection method to achieve a target average network firing rate
F
Setup of the Balanced Networks
The coupling strengths in our networks were scaled with 1/vK, such that the magnitudes of the input current fluctuations had the same order of magnitude in all studied networks. Assuming that inputs from different presynaptic neurons exhibit weak correlations, which is justified in balanced network 13 88 the collective input spike train received by neuron can be effectively modeled as Poisson process with rate denoted as Ω2 = EjEpo(()) Vj = K = Ω. Here, V represents the average firing rate of the network, and K denotes the average number of presynaptic neurons. Assuming the compound input spike train follows a Poisson process, the autocorrelation function of the input current can be expressed as follows:
C(T) (δI(t)δI(t+ + T))e
(14)
) S) -s)ds
Jo
(15)
VR
J3 (57()) K
(16)
((P)(T)
(17)
where 8 denotes the Dirac delta function. In the diffusion approximation, characterized by high Poisson rate and weak coupling, fluctuations in the input currents can be described as Gaussian white noise with magnitude given by
=-
(18)
Note that due the scaling of the coupling strengths J Jo with the square root of the number of synapses K the magnitude of the fluctuations 0² is independent of the number of synapses. Therefore, the input fluctuations do neither vanish nor diverge in the thermodynamic limit and the balanced state in sparse networks emerges robustly [86] 87).
VK
The existence of fixed point of the population firing rate in the balanced state for large follows from the equation of the network-averaged mean current:
I V - (oD)
(19)
In the limit of large K, self-consistency requires the balance between excitation and inhibition, resulting in neurons firing at their maximum rate. The resulting strong inhibition would break the inequality, JoB) < O, the mean current / would diverge to oo, causing the neurons to be silent. Again, the lack of inhibition breaks the inequality. The self-consistency condition in the large K limit can only be satisfied when
specifically To Joh.. If lim (To JoD) > 0, the mean current I would diverge
leading to contradiction. On the other hand, if limm
Lo Jop O
which ensures that the excitatory external drive and the mean recurrent inhibitory current cancel each other. Note that since to Jop O(1/VK /K) the network mean current has a finite value in the large K-limit. The average population firing rate, expressed in units of the membrane time constant Tm can be approximated
Io J0 +
(Vx))
1
(20)
Jo
This approximation generally becomes exact for large K. For excitatoryiinhiiioory mixed network. an analogous self-consistency argument can be made which results in set of inequalities that must be fulfilled achieve balanced state [87]
G Code Availability
An example implementation of networks of LIF, QIF and EIF neurons using SparseProp in Julia, in the balanced state, is available at http://ghthhao.ooonnnniiccegeear
20


Reviews:
[{"review_id": "YERApixJW9", "reviewer": "Reviewer_XBiB", "summary": "This paper proposes SparseProp, an event-based algorithm for efficient simulation and training of sparse Spiking Neural Networks (SNNs). By leveraging the sparsity of the network, SparseProp reduces the computational costs of operations from O(N) to O(log(N)) per network spike. This method avoids iterating through all neurons with each spike and adopts efficient state updates. Demonstrations include a simulation of a sparse SNN with one million LIF neurons, achieving speeds surpassing prior models by over four orders of magnitude. This innovation offers the potential for more sophisticated brain-inspired models.", "questions": "1. Can the proposed simulator scale in a paralleled system? Moreover, how well is the scalability?\n\n2. Have the authors considered using the simulator for layered network architecture? What will the difference be compared with the paper's recurrent architecture?\n\n3. What will be the performance comparison regarding algorithm accuracy and computation cost when comparing an SNN using the event-based simulator and an SNN using a step-based approach running in parallel computing frameworks?", "limitations": "The authors have listed the limitations in the paper. The reviewer agrees with the author that the major limitation of the paper is the need for actual workload and comparison.", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The proposed method decreases the computational cost of conventional event-based SNN simulation from O(n) to O(log(n)) per network spike. The cost decrease mainly comes from using efficient data structures for spike searching and change of reference frame for state evolution.\n\n2. The authors demonstrate the potential to use the event-based simulator for training SNN. This can be a handy tool for SNN training combined with the EventDrop algorithm. It could decrease the training cost of the SNN when number of events is low. However, the authors performed no experiments on this matter, making the contribution vague.", "weaknesses": "The main concern of the reviewer is the usefulness of the proposed approach in real-world applications. Thus, the proposed method may better serve as a neuroscience tool than be used in real applications, as advertised by the authors in the introduction. Here is a list of weaknesses of the paper the reviewer summarized:\n\n1. The paper lacks experiments with SNN algorithms for real-world tasks. The experiments done in the paper concentrate on analytical modeling and non-realistic neural network architecture. However, SOTA SNN algorithms have different architectures and may generate different results. The authors need to perform experiments on real tasks to make a convincing statement on the applicability of real applications.\n\n2. Event-based computation typically suffers from parallelization problems and scalability problems. The problem becomes more apparent when the number of spikes increases in the SNN compared to sparse networks. On the other hand, GPUs and typical deep-learning accelerators use vector computation which is very easy to parallel and can generate much lower latency compared with the event-based approach. The authors must discuss and compare the two approaches using actual workload and show under what conditions the event-based approach may have an advantage.\n\n3. It is hard to say that the proposed method will have an advantage compared with the significant amount of SOTA SNN algorithms that use the Euler method to simulate the activity of spiking neurons. The step-based simulation method demonstrates its effectiveness in many complex tasks and efficiency in many neuromorphic hardware. Therefore, the significance of the proposed approach can be limited.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the thoughtful comments, which substantially improved the manuscript. As requested, we are working towards evaluating SparseProp on real-world scenarios, and we added clarifications on parallelization and comparison to SOTA SNN algorithms that use the Euler method to simulate the activity of spiking neurons.\n\nHere are our more detailed responses to the questions:\n\n>   Can the proposed simulator scale in a paralleled system? Moreover, how well is the scalability?\n\nGreat question. Indeed, parallelization of binary heaps and priority queues is possible [1] but goes beyond the scope of this publication. We would like to emphasize that while parallelization might in an ideal scenario give a linear speedup with the number of threads or processors involved, parallelization does not generally reduce the computational cost in terms of flops per spike or energy cost (e.g., Joule per spike). In contrast, SparseProp mirrors the speedup of massive parallelization (speedup of 10^4 for 10^6 neurons) with much lower computational and energy costs. We leave the proposed parallel heap implementation as a promising future research avenue.\n\n[1] J.-R. Sack and T. Strothotte, An Algorithm for Merging Meaps, Acta Informatica 22, 171 (1985).\n\n> Have the authors considered using the simulator for layered network architecture? What will the difference be compared with the paper's recurrent architecture?\n\nAnother great suggestions. Layered recurrent network architectures can also be simulated using our implementation. We added another Figure and supplementary code to demonstrate this. We chose a 3-layered recurrent network. An extension of SparseProp to feedforward networks is also possible, but this publication is focussing on spiking recurrent networks. We clarified this in the updated manuscript. \n\n> What will be the performance comparison regarding algorithm accuracy and computation cost when comparing an SNN using the event-based simulator and an SNN using a step-based approach running in parallel computing frameworks?\n\nGreat suggestion, we are in the process of benchmarking against NEST, which is a commonly used spiking network simulator."}, {"review_id": "JaMb3sw4Ti", "reviewer": "Reviewer_gRuk", "summary": "This paper presents an efficient event-based algorithm named *SparseProp* for both inference and learning of sparse SNNs. It gets rid of the discretized simulation of ODEs in neuronal dynamics and utilizes phase representation for those neuron models with closed-form solutions between spike times. A priority queue is employed for recording and evolving the states of related neurons when a spike is triggered, which reduces the amortized time complexity to O(K log N) per spike.", "questions": "1. *SparseProp* seems only works for some specific cases, such as the spontaneous firing neuron model. Is it possible to generalize *SparseProp* to a broader variety of SNNs?\n2. Does every operation in *SparseProp* have a well-defined gradient for backpropagation? - The reset of phase for the firing neurons, the phase changes of postsynaptic neurons, and the reset of the global phase (Lines 100-101), to name a few.", "limitations": "See the weaknesses part", "rating": 4, "confidence": 3, "soundness": 1, "presentation": 3, "contribution": 3, "strengths": "This work manages to perform a lazy update of neuron states. Phase representation rather than membrane potential is a brilliant idea since all neurons now share an identical magnitude of update in internal states, which makes it feasible to update $\\phi^{reset},\\phi^{th}$ alone as an alternative to the O(N) polling all neurons' potentials. The priority queue ensures the incoming spike always happens in the neuron at the top of the heap. The experimental results also show significantly improved simulation efficiency. Although there exists earlier work [1] applying priority queue to the simulation of event-based SNNs, the algorithm presented here for sparse SNNs is definitely not trivial.\n\n[1] Eduardo Ros, et al. \"Event-driven simulation scheme for spiking neural networks using lookup tables to characterize neuronal dynamics.\" Neural computation, 2006", "weaknesses": "There are still some major concerns to be resolved\n\n1. The whole study seems to assume a positive external input, i.e. $I^{ext}>0$ in Line 155/544 and a special setting $V_{th}=0, V_{reset}=-1$ in Line 134. In such a case, we expect a spontaneous spike dynamic for those integrate-and-fire models. We have a well-defined unperturbed interspike interval $T^{free}$ and can deduce the phase $\\phi$ accordingly. However, a popular setting in many recent SNN studies, such as ref.6 and ref.17, is $V_{th}>0, V_{reset}=0$, and $I^{ext}=0$. Under this assumption, a neuron won't fire if the incoming excitatory spike is not strong enough. There will be no spontaneous spike and thus no unperturbed interspike interval $T^{free}$ or valid phase representation, and the whole thing seems to fall back to polling of neuron states. In fact, any $I^{ext}\\in[V_{reset}, V_{th}]$ lead to an invalid $T^{free}$ for the LIF model as Eq.7 shows. All in all, I don't think we can assume those parameters without losing generality as Lines 133-134 reads.\n2. For an excitatory-dominated network, it might be possible that two neurons become postsynaptic neurons of each other. *SparseProp* may find that any firing of one will trigger a spike of the other, and it implies an endless loop of updating the priority queue.\n3. The training algorithm coupled with *SparseProp* is more like a draft since no experimental results are really attached here. There might be many details to be considered when designing a backpropagation algorithm for *SparseProp*.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the detailed and insightful feedback and request for clarifications, which substantially improved the publications. \nHere are our more detailed responses to the questions:\n\n> The whole study seems to assume a positive external input, i.e. in Line 155/544 and a special setting in Line 134. In such a case, we expect a spontaneous spike dynamic for those integrate-and-fire models. We have a well-defined unperturbed interspike interval and can deduce the phase accordingly. However, a popular setting in many recent SNN studies, such as ref.6 and ref.17, is , and . Under this assumption, a neuron won't fire if the incoming excitatory spike is not strong enough. There will be no spontaneous spike and thus no unperturbed interspike interval or valid phase representation, and the whole thing seems to fall back to polling of neuron states. In fact, any lead to an invalid for the LIF model as Eq.7 shows. All in all, I don't think we can assume those parameters without losing generality as Lines 133-134 reads.\n\nThanks for allowing us to clarify this. It is accurate that the phase representation of QIF and the pseudo-phase representation of LIF require positive external input currents. **However, our proposed SparseProp algorithm also works with negative input currents.** In that case, the binary heap contains the time to the next spike for each neuron. If there is no recurrent neuron in the network that has a finite next spike time because of negative input currents, one might need input spike  (e.g., Poisson), which can also be handled by the SparseProp. We demonstrate this here for a network of quadratic integrate-and-fire (QIF) neurons with excitatory Poisson input. In that case, one can still solve the single neuron dynamics V(t) analytically between spikes. For dimensionless QIF neuron dynamics\n$\\frac{dV}{dt}=V^2-I_{ext}$ \nwith initial condition $V(0)=V_0$, this would be \n\n$V(t) = -\\sqrt{I_{ext}} \\tanh(\\sqrt{I_{ext}} t - \\tanh^{-1}\\left(V_0/\\sqrt{I_{ext}}\\right))$. \n\nAs soon as an external (or recurrent) input spike pushes one of the recurrent neurons beyond the spike threshold, this neuron has a finite time to the next spike, which leads to the usual treatment in the binary heap.  The time to the next spike for the quadratic integrate-and-fire neuron with $V_0>\\sqrt{I_{ext}}$ is:\n\n$t_{spike} =   \\frac{π - 2\\tanh^{-1}\\left(\\frac{V_0}{\\sqrt{I_{ext}}}\\right)}{2\\sqrt{I_{ext}}}$\n\nWe are working on adding a SparseProp implementation with Poisson input spike trains to the code repository as a demonstration.\n\nNegative external input currents for neuron models that lack an analytical solution $V(t)$ can be handled, as in the example of the exponential integrate-and-fire neuron: One can use Chebyshev polynomials for efficient state update and calculation of the next spike time in SparseProp.\n\n> For an excitatory-dominated network, it might be possible that two neurons become postsynaptic neurons of each other. SparseProp may find that any firing of one will trigger a spike of the other, and it implies an endless loop of updating the priority queue.\n\nIndeed, if two excitatory neurons are postsynaptic to each other, in certain cases a positive feedback loop can occur, which is a problem in any implementation as the firing rate of the two neurons could diverge. Such a scenario can also occur in forward Euler, and glancing over the issue by keeping a fixed integration time step $\\Delta t$ despite the divergence will yield mathematically incorrect results. \nInstead, a solution (that can solve such an issue) would be to either have an inhibition-dominated network and/or a refractory period or synaptic delay such that an excitatory neuron can not instantaneously spike again. Finally, we want to stress that a divergence of the firing rate can be avoided in a network of Quadratic integrate-and-fire neurons, as QIF neurons still take a finite time to spike even after being pushed beyond their firing threshold.\n\n> The training algorithm coupled with SparseProp is more like a draft since no experimental results are really attached here. There might be many details to be considered when designing a backpropagation algorithm for SparseProp.\n\nWe agree with the reviewer that training spiking networks with SparseProp on tasks is an important future direction. What we propose here is to use SparseProp to boost the performance of the already existing event-based EventProp, thus no extra backprop algorithm for SparseProp needs to be designed [1, 2]. However, as pointed out in the limitation paragraph (line 216), a fundamental limitation of directly using SparseProp for training is the apparent incompatibility of event-based simulations with surrogate gradient techniques [3, 4]. While one could introduce 'ghost spikes' in the event-based simulation to emulate surrogate gradients when neurons hit a lower threshold. However, it remains unclear how to preserve the favorable computational scaling of SparseProp in this case. We further clarified this in the training paragraph.\n\n[1] T. C. Wunderlich and C. Pehle, Event-Based Backpropagation Can Compute Exact Gradients for Spiking Neural Networks, Sci Rep 11, 12829 (2021).\n\n[2] T. Nowotny, J. P. Turner, and J. C. Knight, Loss Shaping Enhances Exact Gradient Learning with EventProp in Spiking Neural Networks, 2022.\n\n[3] E. O. Neftci, H. Mostafa, and F. Zenke, Surrogate Gradient Learning in Spiking Neural Networks: Bringing the Power of Gradient-Based Optimization to Spiking Neural Networks, IEEE Signal Processing Magazine 36, 51 (2019).\n\n[4] F. Zenke and S. Ganguli, SuperSpike: Supervised Learning in Multilayer Spiking Neural Networks, Neural Computation 30, 1514 (2018).\n\n\n\n\n"}, {"review_id": "Lnkc1GDJ8u", "reviewer": "Reviewer_h5Gz", "summary": "Simulation of RSNN (recurrent spiking neural network) is a difficult task, whose limits prevent large scale (say 1e.9 neurons) simulations. Often, in practice, one neuron is connected to only a few others, the sparsity. The present paper leverages this sparsity to upscale the simulations limits of RSNN. Sparsity is often leveraged to increase performance and scalibity but SparseProp (algo. presented in this paper) seems to be a newcomer for event-based simulations of RSNN.", "questions": "-What is the interrogation point, line 190 ? Probably a reference to correct\n-Which type of network is used in Figure3 ? Is the scaling the same for heterogenous networks ?", "limitations": "-limitations raised by the the weaknesses / questions stated earlier.\n-limitations raised by the authors themselves in section 8", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "-The authors present a new and efficient algorithm for event-based simulations of RSNN. In its structure, the algorithm is fairly similar to classical ones, \nwhile taking advantage of the sparsity structure, it presents a very good ratio (high) efficiency / complexity. Practically, it mainly forces the practitionner to have a better data structure.\n-The exposition and background are in general clear and the algorithm seems natural to derive in the history-flaw of the paper (great!). Algorithms for non-sparse and sparse simulations\nare cleary written and can be easily compared by the reader.\n-Simulations are given, for SparseProp and classical algo, giving a good sense on how performant their method is (performant!).\n-The flow of the paper Abstract / General -> Particular / Practical, is smooth.\n-Examples on classical RSNN are given, which are good toy examples for the reader to play with.\n-Simulations are available in Julia.", "weaknesses": "-General comment on citations: the authors could give more references, for example:\n\t-line 59: several classical papers showing the origin of the differential equation, would be appreciated.\n\t-figure 1: a reference paper, explaining the role and subtilities of the phase representation, would be appreciated (especially since its of crucial importance for this paper).\n-The claim of efficient simulation is for general univariate RSNN but the paper nearly entirely focuses on homogeneous ones. I understand that this can be generalized, \nas explained in Section 6, but no numerical results (as in Figure 3) are given and the algorithm is only explained (Section 2.1) for homogeneous networks. For pedagological purpose, \nits great, but it leaves with a taste of unfinished. It feels like the most interesting and difficult part is \"left for the reader\". I would suggest claiming a weaker result or \ngiving more material for the understanding of heterogeneous networks (eventually in the appendix). This is the main reason why Soundness grade is only at 2.\n-Since the pulse-coupled phase oscillators (and constant phase velocity) is an important and simplifying assumption, the paper lacks a discussion on it or a reference that discuss it (is it a strong assumptions: yes, no, why - how is it linked to Eq. 1...).\n-Section 7 is simulation-less so the reader can not know if the proposed algorithm is efficient. Moreover it is mentionned \"training of RSNN\" in the title, which seems a bit of an overstatement since there is only a 20 line, simulation-less, section about it. I know that this is because the authors of EventProp did not make their algo available but this also should not be a problem for the reader. I would suggest removing the \"training\" in the title and keeping this section either in an appendix or as a \"bonus\" section. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time and helping us improve our work. We will integrate these ideas into our manuscript with care. Here are our detailed responses:\n\n> General comment on citations: the authors could give more references, for example: -line 59: several classical papers showing the origin of the differential equation, would be appreciated.\n\nFixed in the updated manuscript.\n\n> figure 1: a reference paper, explaining the role and subtilities of the phase representation, would be appreciated (especially since its of crucial importance for this paper). \n\nWe added several references to the updated manuscript. While the phase representation is being used in the examples in the main manuscript, we want to stress that SparseProp also can be used in neurons **without** phase representation (e.g., exponential integrate-and-fire neurons) and in neurons, that receive negative external input.\n\n> The claim of efficient simulation is for general univariate RSNN but the paper nearly entirely focuses on homogeneous ones. I understand that this can be generalized, as explained in Section 6, but no numerical results (as in Figure 3) are given and the algorithm is only explained (Section 2.1) for homogeneous networks. For pedagological purpose, its great, but it leaves with a taste of unfinished. It feels like the most interesting and difficult part is \"left for the reader\". I would suggest claiming a weaker result or giving more material for the understanding of heterogeneous networks (eventually in the appendix). This is the main reason why Soundness grade is only at 2. \n\nWe appreciate the suggestions for improvements. However, we want to point out that the supplement C already contains pseudocode for a heterogeneous SparseProp and mathematical details for their implementation with quadratic integrate-and-fire models. Furthermore, we provide code for the heterogeneous SparseProp implementation. Thus, we would like to politely disagree that \"details are left to the reader\". We however do agree that from the main manuscript the details of the heterogeneous SparseProp does not become fully clear. We plan to use the additional page for the camera-ready version for a more pedagogical explanation, should the paper be accepted.\n\n> Since the pulse-coupled phase oscillators (and constant phase velocity) is an important and simplifying assumption, the paper lacks a discussion on it or a reference that discuss it (is it a strong assumptions: yes, no, why - how is it linked to Eq. 1...). \n\nWe added references that discuss this to the updated manuscript. Overall, as argued in the paragraph on limitations, one core limitation of event-based simulations is that one needs a solution of $V(t)$ -either in closed form or as a lookup-table. That means that incorporating stochastic input without an analytical solution of $V_i(t)$ is in the current implementation not possible. We would like to stress that we introduce several extensions in the supplement and in the provided online code that go beyond the pedagogically useful but restrictive LIF and QIF implementation of the main text:\n* SparseProp for heterogeneous networks with a time-based binary (which would for example allow different constant external input current and different $\\tau_M$ for each neuron. \n* SparseProp for neuron models that lack an analytical closed-form expression for the next spike time, like the exponential integrate-and-fire model.\n* SparseProp with Poisson input. For Poisson input, the speedup of SparseProp becomes particularly relevant, as in the conventional implementation every Poisson input spike would be of computational complexity $O(N)$.\n\n\n> Section 7 is simulation-less so the reader can not know if the proposed algorithm is efficient. Moreover it is mentionned \"training of RSNN\" in the title, which seems a bit of an overstatement since there is only a 20 line, simulation-less, section about it. I know that this is because the authors of EventProp did not make their algo available but this also should not be a problem for the reader. I would suggest removing the \"training\" in the title and keeping this section either in an appendix or as a \"bonus\" section.\n\nWe agree with the reviewer that the title might raise wrong expectations, and we will follow the reviewer advice and remove the “training” from the title.\n\nQuestions:\n\n> What is the interrogation point, line 190 ? Probably a reference to correct\n\nYes, we fixed this typo. The missing references were:\n[1] D. Liberzon, Calculus of Variations and Optimal Control Theory: A Concise Introduction, in Calculus of Variations and Optimal Control Theory (Princeton University Press, 2011).\n[2] L. S. Pontryagin, Mathematical Theory of Optimal Processes: The Mathematical Theory of Optimal Processes, 1st edition (Routledge, New York, 1987).\n\n\n\n > Which type of network is used in Figure3 ? \n\nFigure 3 used a recurrent network of quadratic-integrate-and fire neurons. We added this now to the caption.\n\n\n > Is the scaling the same for heterogenous networks ?\n\nYes, as demonstrated in the supplement and in the available online code, the favorable scaling of $O(log(N))$ of SparseProp compared to O(N) for the conventional implementation is maintained for heterogeneous networks."}, {"review_id": "BkuNgZ7Sab", "reviewer": "Reviewer_MMBd", "summary": "This paper introduces an efficient Spike-based recurrent network learning simulation method based on an event-driven implementation.\n\nThe method is based on two main contributions. The first consists in a change of temporal reference frame that uses the phase of a neuron rather than its absolute time, the second contribution consists in finding an efficient data representation corresponding to the parsimonious structure of neural connectivity. This method is applied to different neuron models, from a linear LIF model to more complex models such as the quadratic LIF neuron.", "questions": "I am curious to know if the change in representation from absolute time to phase affects the precise latency of neuron firing? Additionally, it appears that biological neural networks include transmission delays and spikes between neurons. Is your implementation suitable for this scenario?\n", "limitations": "Since the computational results seems so promising, a limit of this paper is that it is not applied to standard benchmarks and notably to datasets which include a rich dynamics.\n", "rating": 9, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "This paper clearly presents an original method for simulating and training a spiking neural network. The motivation for the paper is very well set out, and the methods are clearly explained both in terms of the definition of the mathematical equations and the path that leads to its numerical implementation. The results look very promising, and the preliminary results set out in the appendix are very promising indeed. In particular, the experimental results show that this method can be used to implement very large networks on computational architectures of reasonable size, which is a game changer for the SNN community.\n\nMinor:\ncaption of Fig 4 : repetition of « higher densities » ", "weaknesses": "Although Figure 3 already presents a comparison between the proposed method and a conventional method, it would have been desirable to have a comparison between these two methods in terms of the quality of the obtained results. It would have also been beneficial to isolate the improvements achieved by the two contributions independently (the first contribution being the change in temporal reference, and the second being the sparse representation of data). Some improvements in the presentation can be made... For example, the \"peek\" function is not defined.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their insightful and encouraging assessment and agree that SparseProp could be a game changer for the spiking network community, especially, for very large, sparse networks.\n\nHere are our more detailed responses:\n\n> Minor: caption of Fig 4 : repetition of « higher densities »\n\nFixed in the revised version of the manuscript.\n\n> Although Figure 3 already presents a comparison between the proposed method and a conventional method, it would have been desirable to have a comparison between these two methods in terms of the quality of the obtained results. \n\nGreat suggestion!  We are  adding code fore a “work-precision diagram” to the anonymous code. This compares the computational cost measured in terms of time per network spike against the accuracy or “precision” of SparseProp, allowing for a comparative assessment of different methods in terms of their efficiency and accuracy.\n\n> It would have also been beneficial to isolate the improvements achieved by the two contributions independently (the first contribution being the change in temporal reference, and the second being the sparse representation of data). \n\nWe were considering that, but as both ingredients are necessary for going from a computational complexity of O(N) to O(log(N)), we decided against that. However, we are working on adding this comparison to the supplement.\n\n> Some improvements in the presentation can be made... For example, the \"peek\" function is not defined.\n\nGood catch, we added a definition now. The “peek” function in a binary heap returns the element with the highest priority without removing it from the heap. In the case of a min-heap, which is used for the phase representation, it retrieves the smallest element (next spike time), while in a max-heap (used for heterogeneous networks), it retrieves the largest element (maximum phase). This operation typically takes constant time, O(1), since the element is at the root of the heap. However, it doesn't modify the heap structure in any way.\n\n> I am curious to know if the change in representation from absolute time to phase affects the precise latency of neuron firing?\n\nGreat question. On paper, the expression is exact and thus does not affect the latency. Numerically, as over the course of many spikes, the global phase offset $\\Delta \\phi$ is incrementally growing, for very large $\\Delta \\phi$, catastrophic cancellation can occur (for example, when $\\Delta \\phi=10.0^6$, and it is subtracted from $d\\phi$ which is very small. In this case, the problem one runs into is that the next bigger floating point number for double floats next to 10.0e6 has a distance of $2.0^{-33}\\approx$1.1641532182693481e-10.  This source of numerical error can be reduced by resetting the global phase shift and all individual neuron's phases when it exceeds some threshold to avoid numerical errors resulting from subtractive cancellation due to floating-point arithmetic. However, this reset has only to occur very infrequently, on the order of every $c N$ steps.  $c$ is a big prefactor that depends on the acceptable floating-point error, e.g., for an error $2.0^{-33}$ it would be 10.0^6, which can be easily checked in Julia by the command:\n\n```julia\nlog2(nextfloat(1.0e6)-1.0e6)\n```\n\n> Additionally, it appears that biological neural networks include transmission delays and spikes between neurons. Is your implementation suitable for this scenario?\n\nYes, thanks for the opportunity to clarify this. Indeed, SparseProp can easily be extended to transmission delays in case of the heterogeneous implementation. In that case, any network spike would not trigger an instantaneous update of the phase response curve of postsynaptic neurons. Instead, the current time plus the synaptic delay would be added as an event to the binary heap and processed at that later time. \n\n> Since the computational results seems so promising, a limit of this paper is that it is not applied to standard benchmarks and notably to datasets which include a rich dynamics.\n\nWe agree that it is an important next step to further benchmark the algorithm. Which benchmarks would the reviewer suggest? We are working on benchmarking the algorithm itself against NEST, but other suggestions would be very welcome."}]

Author rebuttal:
We deeply appreciate the feedback and constructive critiques provided on our submission. Here, we provide a consolidated response to address the concerns and suggestions raised.

1. **Reference and Citations:** In response to the recommendation for added references, we've updated the manuscript to include pivotal papers related to the differential equation origins, phase representation, and the role of pulse-coupled phase oscillators. We believe these additions fortify the foundational knowledge upon which our work stands.

2. **Clarity on Homogeneous vs. Heterogeneous Networks:** We acknowledge the pointed out discrepancy between the focus on homogeneous networks and the broader claim of the algorithm's applicability. We would like to emphasize that the supplementary material C provides pseudocode and mathematical details for implementing SparseProp in heterogeneous contexts. The paper's digital resources also offer code for the heterogeneous implementation. To further enhance clarity, we've allocated additional space in the camera-ready version to elucidate the workings of the heterogeneous SparseProp.

3. **Soundness of the Algorithm and Result Presentation:** We concur with the reviewer's observation regarding the seeming lack of simulation in Section 7. To address this, we've refined the title to better align with content expectations. The typos, specifically the interrogation point in line 190, have been rectified, and Figure 3 now clearly states the type of network in its caption. Moreover, we've added comparative elements showcasing the computational cost versus the accuracy of SparseProp, highlighting its efficiency.

4. **Specific Technical Queries:** The "peek" function, which was previously ambiguous, has now been defined explicitly. The question on the influence of representation change on neuron firing latency is pertinent. The paper emphasizes that, theoretically, the expression remains unaffected, but numerical considerations are taken into account to minimize errors. Furthermore, our methodology is flexible enough to incorporate transmission delays, offering more comprehensive applicability. Regarding the consideration of real-world benchmarks, we appreciate suggestions on this front and are exploring benchmarking against established platforms like NEST.

5. **Additional Enhancements:** In light of the feedback received, we added an extra page elucidating nontrivial network architectures whose dynamics go beyond asynchronous irregular activity. We've also incorporated a comparative study of the numerical phase response curve and solutions based on Chebyshev polynomials. A new figure has been introduced to visualize the error in individual spike times of SparseProp in a Leaky-integrate and fire neuron network as requested by one reviewer. 

We believe these modifications and clarifications have significantly strengthened our submission, aligning it more closely with NeurIPS standards. We once again extend our gratitude for the insightful feedback and hope our revisions will cast our work in a more favorable light.

Meta-review: This paper proposes a novel method to substantially scale up recurrent spiking neural networks by exploiting sparsity in the connection matrix. It elicited a very high level of enthusiasm from two reviewers, who praised its clarity, originality, and thoroughness of exposition, and two more negative appraisals from reviewers, who raised concerns about novelty and originality of SparseProp and its usefulness for real-world applications.  Having read the detailed reviews and rebuttals, it seems to me that the enthusiasm of the positive reviewers outweighed the weaknesses pointed out by the more critical reviewers, and I believe the paper should be accepted to this year's NeurIPS.  Please make a thorough attempt to address the reviewer comments (particularly those of reviewer gRuK about clarity and soundness) in the final manuscript. Congratulations!
Decision: Accept (poster)