Event-based backpropagation can compute exact gradients for spiking neural networks

Spiking neural networks combine analog computation with event-based communication using discrete spikes. While the impressive advances of deep learning are enabled by training non-spiking artificial neural networks using the backpropagation algorithm, applying this algorithm to spiking networks was previously hindered by the existence of discrete spike events and discontinuities. For the first time, this work derives the backpropagation algorithm for a continuous-time spiking neural network and a general loss function by applying the adjoint method together with the proper partial derivative jumps, allowing for backpropagation through discrete spike events without approximations. This algorithm, EventProp, backpropagates errors at spike times in order to compute the exact gradient in an event-based, temporally and spatially sparse fashion. We use gradients computed via EventProp to train networks on the Yin-Yang and MNIST datasets using either a spike time or voltage based loss function and report competitive performance. Our work supports the rigorous study of gradient-based learning algorithms in spiking neural networks and provides insights toward their implementation in novel brain-inspired hardware.

We derive the precise analogue to backpropagation for spiking neural networks by applying the adjoint method together with the jump conditions for partial derivatives at state discontinuities, yielding exact gradients with respect to loss functions based on membrane potentials or spike times.A, B: Dynamical systems with parameter-dependent discontinuous state transitions typically have discontinuous partial derivatives of state variables with respect to system parameters (1), as is the case for the two examples shown here.
Both examples model dynamics occurring on short timescales, namely inelastic reflection and the neuronal spike mechanism, using an instantaneous state transition.We denote quantities evaluated before and after a given transition by − and +.In A, a bouncing ball starts at height y 0 > 0 and is described by ÿ = −g with gravitational acceleration g.It is inelastically reflected as ẏ+ = −0.8ẏ− as soon as y − = 0 holds, causing the partial derivative with respect to y 0 to jump as ∂y + ∂y0 = −0.8∂y − ∂y0 (see section 4.1).In B, a leaky integrate-and-fire neuron described by the system given in table 1 with initial conditions I(0) = w, V (0) = 0 resets its membrane potential as V + = 0 when V − = ϑ holds, causing the partial derivative to jump as ∂w (see section 4.3).C: Applying the adjoint method with partial derivative jumps to a network of leaky integrate-and-fire neurons (table 1) yields the adjoint system (table 2) that backpropagates errors in time.EventProp is an algorithm (algorithm 1) returning the gradient of a loss function with respect to synaptic weights by computing this adjoint system.The forward pass computes the state variables V (t), I(t) and stores spike times t post and each firing neuron's synaptic current.EventProp then performs the backward pass by computing the adjoint system backwards in time using event-based error backpropagation and gradient accumulation: each time a spike was transferred across a given synaptic weight in the forward pass, EventProp backpropagates the error signal represented by the adjoint variables λ V (t post ), λ I (t post ) of the post-synaptic (target) neuron and updates the corresponding component of the gradient by accumulating λ I (t post ), finally yielding sums as given in the figure.

Introduction
How can we train spiking neural networks to achieve brain-like performance in machine learning tasks?The resounding success and pervasive use of the backpropagation algorithm in deep learning suggests an analogous approach.This algorithm computes the gradient of the neural network parameters with respect to a loss function that measures the network's performance in a given task.The parameters of the network are iteratively updated using the locally optimal direction given by the gradient.Spiking neural networks have been referred to as the third generation of neural networks (2), superseding artificial neural networks as commonly used in deep learning and hold the promise for efficient and robust processing of event-based spatio-temporal data as found in biological systems.However, spiking models are not widely used in machine learning applications.At the same time, the development of spiking neuromorphic hardware receives increasing attention (3) and learning in spiking neural networks is an active research subject, with a wide variety of proposed algorithms.A notorious issue in spiking neurons is the hard spiking threshold that does not permit a straight-forward application of differential calculus to compute gradients.Although exact gradients have been derived for special cases, this issue is commonly side-stepped by using smoothed or stochastic neuron models or by replacing the hard threshold function using a surrogate function, leading to the computation of surrogate gradients (4).
In contrast, this work provides an algorithm, EventProp, to compute the exact gradient for an arbitrary loss function defined using the state variables (spike times and membrane potentials) of a general recurrent spiking neural network composed of leaky integrate-and-fire neurons with hard thresholds.Since feed-forward architectures correspond to recurrent neural networks with block weight matrices and convolutions can be represented as sparse linear transformations, deep feed-forward networks and convolutional networks are included as special cases.

Partial Derivatives in Discontinuous Dynamical Systems
The leaky integrate-and-fire neuron model describes a hybrid dynamical system that combines continuous dynamics between spikes with discontinuous state variable transitions at spike times.The computation of partial derivatives for hybrid dynamical systems is an established topic in optimal control theory (1; 5).In hybrid systems, the time-dependent partial derivative ∂x ∂p (t) of a state variable x with respect to a parameter p generally experiences jumps at the points of discontinuity (see fig. 1 A, B).The relation between the partial derivatives before and after a given discontinuity was first studied in the 1960s (6; 7).A more general theoretical framework was developed thirty years later (8), providing existence and uniqueness theorems for the partial derivative trajectories ∂x ∂p (t) of hybrid systems.Discontinuous state transitions in hybrid systems occur when a transition condition is fulfilled (e.g., a bouncing ball hits the floor or a neuron reaches its spiking threshold).The existence of well-defined partial derivative jumps at the state transition times depends on the local applicability of the implicit function theorem to the transition condition, requiring that the event time depends on the parameters in a differentiable fashion.In the case considered here, a spiking neural network composed of leaky integrate-and-fire neurons that is parameterized by synaptic weights, this is fulfilled up to the null set in weight space that contains the locally defined hypersurfaces where spikes are added or removed.At these critical points, the derivative of the time of the (dis-)appearing spike with respect to a given active synaptic weight diverges.This implies that both the spike times and an integral of a smooth loss function over the membrane potential are differentiable almost everywhere, up to the null set of critical points in weight space.

Backpropagation of Errors in Discontinuous Dynamical Systems
Having established the jumps of partial derivatives in the leaky integrate-and-fire neuron model, the relevant question is how to compute the gradient of a loss function for spiking neural networks, preferably with the computational efficiency afforded by the backpropagation algorithm and retaining any potential advantages of event-based communication.Backpropagation in discrete-time artificial neural networks can be derived as a special case of the adjoint method (9), with the adjoint variables (Lagrange multipliers) λ t at each time step t corresponding to the intermediate variables computed in the backpropagation algorithm.Applying the adjoint method to continuous-time dynamical systems yields time-dependent adjoint variables λ(t) (see section 4.2) and their computation in reverse time is analogous to the backpropagation of errors in discrete-time artificial neural networks.The adjoint method can be applied to hybrid systems by using the proper partial derivative jumps that generally cause jumps in the adjoint variables (10).

EventProp: Event-Based Backpropagation of Errors
We combine the partial derivative jumps of the leaky integrate-and-fire neuron with the adjoint method in order to derive the EventProp algorithm (algorithm 1) that is the analogue to backpropagation for spiking neural networks (fig. 1 C).Since EventProp backpropagates errors at spike times, the algorithm computes gradients using an event-based communication scheme and is amenable to neuromorphic implementation.By requiring the storage of state variables only at spike times, it provides favorable memory requirements compared to approaches that require the full forward state trajectory to be retained for the backward pass.
For example, surrogate gradient approaches operating on a discrete time grid typically require storing state variables at every time step for the backward pass (however, an effort to compute surrogate gradients in a more sparse fashion has been made 11).More generally, the fact that backpropagation in discrete-time artificial neural networks requires storing activations at every time step causes a memory bottleneck and is a major concern in training very deep architectures (e.g., 12; 13; 14).
EventProp does not prescribe a specific numerical scheme to compute state variables and spike times but since the backward pass corresponds to the computation of a spiking network with pre-determined spike times, the computational complexity of the backward pass generally corresponds to that of the forward pass.While surrogate gradient approaches on a discrete time grid typically require the calculation of dense matrix-vector products at every time step in the backward pass (all neurons backpropagate error signals at every time step), EventProp only requires computing vector-vector products at spike events (only the firing neuron receives backpropagated errors at a given spike time).In this way, EventProp leverages the sparseness of spike-based communication for both the forward and backward pass.
We demonstrate the training of spiking neural networks with a single hidden layer using EventProp and the Yin-Yang and MNIST datasets, resulting in competitive classification performance.

Previous Work
We refer the reader to the following review articles for a comprehensive survey of gradient-based approaches to learning in spiking neural networks: (15) and ( 16) discuss learning in deep spiking networks, (3) discuss learning along with the history and future of neuromorphic computing and (4) focus on the surrogate gradient approach.Surrogate gradients use smooth activation functions for the purposes of backpropagation and have been used to train spiking networks in a variety of settings (e.g., 17; 18; 19; 20).This approach is typically derived by considering the Euler discretization of a spiking neural network where the Heaviside step function is used to couple neurons across discrete time steps.The non-differentiable Heaviside step function is then replaced by a smooth function in the backward pass.
Apart from surrogate gradients, several publications provide exact gradients for first-spike-time based loss functions and leaky integrate-and-fire neurons: (21) provides the gradient for at most one spike per layer and this result was subsequently generalized to an arbitrary number of spikes as well as recurrent connectivity (22; 23).While these publications provide recursive relations for the gradient that can be implicitly computed using backpropagation, we explicitly provide the dynamical system that implements backpropagation through time and show that it represents an adjoint spiking network which transmits errors at spike times, allowing for an event-based computation of the gradient.In addition, we also consider voltage-dependent loss functions and our methodology can be applied to neuron models without analytic expressions for the post-synaptic potential kernels.
The applicability of methods from optimal control theory (i.e., partial derivative jumps and the adjoint method) to compute exact gradients in hard-threshold spiking neural networks was recognized in a series of publications (24; 25; 26).In contrast to this work, these articles consider a neuron model with a two-sided threshold (including negative threshold crossings), rely on the existence of analytic expressions for the post-synaptic potential kernels, provide specialized algorithms tailored to specific loss functions and consider minimalistic regression tasks.
The chronotron (27) uses a gradient-based learning rule based on the Victor-Purpura metric which enables a single leaky integrate-and-fire neuron to learn a target spike train.Our work, as well as the works mentioned above which derive exact gradients, applies the implicit function theorem to differentiate spike times with respect to synaptic weights.A different approach is to consider ratios of the neuronal time constants where analytic expressions for first spike times can be given and to derive the corresponding gradients, as done in (28; 29; 30; 31).Our work encompasses the contained methods to compute the gradient as special cases.
The seminal Tempotron model uses gradient descent to adjust the sub-threshold voltage maximum in a single neuron (32) and has recently been generalized to the spike threshold surface formalism (33) that uses the exact gradient of the critical thresholds ϑ * k at which a leaky integrate-and-fire neuron transitions from emitting k to k − 1 spikes; computing this gradient is not considered in this work.The adjoint method was recently used to optimize neural ordinary differential equations (34) and neural jump stochastic differential equations (35) as well as to derive the gradient for a smoothed spiking neuron model without reset (36).

Results
We first define the used spiking neuron model and then proceed to state our main results.

Leaky Integrate-and-Fire Neural Network Model
We define a network of N leaky integrate-and-fire neurons with arbitrary (up to self-connections) recurrent connectivity (table 1).We set the leak potential to zero and choose parameter-independent initial conditions.Note that the Spike-Response Model (SRM) (37) with double-exponential or α-shaped PSPs is generally an integral expression of the model given in table 1 with corresponding time constants.

Free Dynamics
Transition Condition Jumps at Transition Table 1: The leaky integrate-and-fire spiking neural network model.Inbetween spikes, the vectors of membrane potentials V and synaptic currents I evolve according to the free dynamics.When some neuron n ∈ [1.
.N ] crosses the threshold ϑ, the transition condition is fulfilled, causing a spike.This leads to a reset of the membrane potential as well as post-synaptic current jumps.W ∈ R N ×N is the weight matrix with zero diagonal and e n ∈ R N is the unit vector with a 1 at index n and 0 at all other indices.We use − and + to denote quantities before and after a given spike.

Gradient via Backpropagation
Consider smooth loss functions l V (V, t), l p (t post ) that depend on the membrane potentials V , time t and the set of post-synaptic spike times t post .The total loss is given by Our main result is that the derivative of the total loss with respect to a specific weight w ji = (W ) ji that connects pre-synaptic neuron i (the firing neuron) to post-synaptic neuron j (the receiving neuron) is given by a sum over the spikes caused by i, where λ I is the adjoint variable (Lagrange multiplier) corresponding to the synaptic current I. Equation ( 2) therefore samples the post-synaptic neuron's adjoint variable (λ I ) j at the spike times caused by neuron i.
After the neuron dynamics given by table 1 have been computed from t = 0 to t = T , the adjoint state variable λ I is computed in reverse time (i.e., from t = T to t = 0) as the solution of the system of adjoint equations defined in table 2. The dynamical system defined by table 2 is the adjoint spiking network to the leaky integrate-and-fire network (table 1) which backpropagates error signals at the spike times t post .
Equation ( 2) and table 2 suggest a simple algorithm, EventProp, to compute the gradient (algorithm 1).Notably, if the loss is voltage-independent (i.e., l V = 0), the backward pass of the algorithm requires only the spike times t post and the synaptic current of the firing neurons at their respective firing times to be retained from the forward pass.The membrane potential at spike times is fixed to the threshold ϑ and therefore implicitly retained; the synaptic current therefore determines the temporal derivative of the membrane potential at the spike time, V − , and needs to be stored for the backward pass.The memory requirement of the algorithm scales as O(S), where S is the number of post-synaptic spikes in the network.A feed-forward architecture corresponds to a block matrix W with each block being a strictly triangular matrix that connects two given layers.In that case, the forward and backward pass can be computed in a layer-wise fashion.
In case of a voltage-dependent loss l V = 0, the algorithm has to store the non-zero components of ∂l V ∂V along the forward trajectory.The loss l V may depend on the voltage at a discrete time t i using the Dirac delta,

Jump at Transition
The adjoint spiking network to table 1 that computes the adjoint variable λ I needed for the gradient (eq.( 2)).The adjoint variables are computed in reverse time (i.e., from t = T to t = 0) with = − d dt denoting the reverse time derivative.(λ − V ) n(k) experiences jumps at the spikes times t post k , where n(k) is the index of the neuron that caused the kth spike.Computing this system amounts to the backpropagation of errors in time.The initial conditions are λ V (T ) = λ I (T ) = 0 and we provide λ − V in terms of λ + V because the computation happens in reverse time.
mem at time t i .Note that in many practical scenarios as found in deep learning, the loss l V depends only on the state of a constant number of neurons, irrespective of network size.If l V depends on the voltage of non-firing readout neurons, we have l + V = l − V and the corresponding term in the jump given in table 2 vanishes.
If l V is either zero or depends only on voltages at discrete points in time, EventProp can be computed in a purely event-based manner.Figure 2 illustrates how EventProp computes the gradient of a spike time based loss function for two leaky integrate-and-fire neurons where one neuron receives Poisson spike trains via 100 synapses and is connected to the other neuron via a single feed-forward weight w.
Algorithm 1 EventProp: Algorithm to compute eq. ( 2).Require: Input spikes, losses l p , l V , parameters W , τ mem , τ syn , initial conditions V (0), I(0) grad ← 0 Compute neuron state trajectory (table 1) from t = 0 to t = T : Forward pass for all spikes k, store spike time t post k and the firing neuron's component of I(t post k ) if l V = 0, also store ∂l V ∂V Compute adjoint state trajectory (table 2) from t = T to t = 0: Backward pass accumulate grad ji ← grad ji − τ syn (λ I ) j for each spike from neuron i to j return grad The lower neuron backpropagates its error signal λ V − λ I at the upper neuron's spike times (indicated by arrows).F, G: Accumulated gradient for one of the 100 input weights of the upper neuron and the weight w connecting the upper and lower neuron.EventProp computes the adjoint variables from t = T to t = 0 and accumulates the gradients by sampling −τ syn λ I when spikes are transmitted across the respective weight.The gradients computed in this way match the gradients computed via central differences (dashed lines) up to a relative deviation of less than 10 −7 .

Simulation Results
We demonstrate learning using EventProp using a custom event-based simulator and the Yin-Yang (38) and MNIST (39) datasets.In both cases, we use a single hidden layer and spike latency encoding of the input data.The Yin-Yang dataset is classified using the time to first spike of a layer of readout neurons while the MNIST dataset is classified using the voltage maxima of a layer of non-firing readout neurons.The simulator computes gradients using EventProp as described in algorithm 1; specifically, it uses an event queue and root-bracketing to compute post-synaptic spike times in the forward pass (using exact integration of the membrane potential, ( 40)) and backpropagates errors by attaching error signals to spikes in the backward pass and using reverse traversal of the event queue.We optimized synaptic weights using the calculated gradients via the Adam optimizer (41), without clipping gradients.
By initializing synaptic weights such that the network started in a non-quiescent state, we found that no explicit regularization of firing rates was needed to obtain the reported results in both cases.Hyperparameters were optimized using Gaussian process optimization (42) and manual tuning using the validation set of the respective dataset.The resulting parameters (see table 4) were then evaluated using the test set.The Yin-Yang dataset ( 38) is a two-dimensional non-linearly separable dataset, with a shallow classifier achieving around 64% accuracy, and it therefore requires a hidden layer and backpropagation of errors for high classification accuracy.Consider that in contrast, the MNIST dataset can be classified using a linear classifier with at least 88% accuracy (39).

Yin-Yang Dataset
Each two-dimensional data point of the dataset (x, y) was transformed into four dimensions as (x, 1−x, y, 1−y) and encoded using spike latencies in the interval [0, t max ] (see fig. 3 D).We added a fixed bias spike at time t bias for a total of five input spikes per data point.The resulting spike patterns were used as input to a two-layer network composed of leaky integrate-and-fire neurons.The output layer consisted of three neurons that each encoded one of the three classes, with each data point being assigned the class of the neuron that fired the earliest spike.
In analogy to (28), we used a cross-entropy loss defined using the first output spike times per neuron, where t post i,k is the first spike time of neuron k for the ith sample, l(i) is the index of the correct label for the ith sample, N batch is the number of samples in a given batch and τ 0 and τ 1 are hyperparameters of the loss function.The first term corresponds to a cross-entropy loss function over the softmax function applied to the negative spike times (we use negative spike times as the class assignment is determined by the smallest spike time) and encourages an increase of the spike time difference between the label neuron and all other neurons.As the first term depends only on the relative spike times, the second term is a regularization term that encourages early spiking of the label neuron.

MNIST Dataset
where V k (t) is the voltage trace of the kth readout neuron, l(i) is the index of the correct label for the ith sample and N batch is the number of samples in a given batch.Note that we can write the maximum voltage as max t V k (t) = V k (t)δ(t − t max )dt with the time of the maximum t max and the Dirac delta δ, allowing us to apply the chain rule to find the jump of λ V k (cf.table 2) at time t max (terms containing the distributional derivative of δ are always zero).
During training, input spikes were dropped with probability p drop in order to avoid overfitting.To obtain a validation set, we extracted and removed 5000 samples from the training set.
Training results are shown in fig. 4.After training, the test accuracy was (97.6 ± 0.1) % (mean and standard deviation over 10 different random seeds).This represents competitive classification performance when compared with previously published results using spiking networks with a single, fully connected hidden layer (table 3).

Discussion
We have derived and provided an algorithm (EventProp) to compute the gradient of a general loss function for a spiking neural network composed of leaky integrate-and-fire neurons.The parameter-dependent spike discontinuities were treated in a well-defined manner using the adjoint method in combination with partial derivative jumps, without approximations or smoothing operations.EventProp uses the resulting adjoint spiking network to backpropagate errors in order to compute the exact gradient.Its forward pass requires computing the spike times of pre-synaptic neurons that transmit spikes to post-synaptic neurons, while the backward pass backpropagates errors at these spike times using the reverse path (i.e., from post-synaptic to pre-synaptic neurons).The rigorous treatment of spike discontinuities in combination with an event-based computation of the exact gradient represent a significant conceptual advance in the study of gradient-based learning methods for spiking neural networks.
An apparent issue with gradient descent based learning in the context of spiking networks is that the magnitude of the gradient diverges at the critical points in parameter space (note the v−1 term in the jump term given in table 2; this term diverges as the membrane potential becomes tangent to the threshold and we have v → 0).Indeed, this is a known issue in the broader context of optimal control of dynamical systems with parameter-dependent state transitions (1; 8).While this divergence can be mitigated using gradient clipping in practice, exact gradients of commonly considered loss functions lead to learning dynamics that are ignorant with respect to these critical points and are therefore unable to selectively recruit additional spikes or dismiss existing spikes.In contrast, surrogate gradient methods continuously transmit errors across neurons and combine these with a non-linear function of the distance of the membrane potential to the threshold.It is therefore plausible that surrogate gradients represent a form of implicit regularization.(4) reports that the surrogate gradient approximates the true gradient in a minimalistic binary classification task while at the same time remaining finite and continuous along an interpolation path in weight space.Hybrid algorithms that combine the exact gradient with explicit regularization techniques could be a direction for future research and provide more principled learning algorithms as compared to ad-hoc replacements of threshold functions.
This work is based on the widely used leaky integrate-and-fire neuron model.Extensions to this model, such as fixed refractory periods, adaptive thresholds or multiple compartments can be treated in an analogous way (47).While the absence of explicit solutions to the resulting differential equations can require the use of sophisticated numerical techniques for event-based simulations, such extensions can significantly enhance the computational capabilities of spiking networks.For example, (18) uses adaptive thresholds to implement LSTM-like memory cells in a recurrent spiking neural network.
Neuromorphic hardware is an increasingly active research subject (e.g., 48; 49; 50; 51; 52; 53; 54; 55; 56; 57; 58) and implementing EventProp on such hardware is a natural consideration.The adjoint dynamics as given in table 2 represent a type of spiking neural network which, instead of spiking dynamically, transmits errors at fixed times t post that are scaled with factors v−1 retained from the forward pass.Therefore, a neuromorphic implementation could store spike times and scaling factors locally at each neuron, where they could be combined with the dynamic error signal (λ V − λ I in table 2) in the backward pass.This requires a possibility to read out neuronal state variables both in the forward and backward pass (membrane potential and synaptic current).The resulting error signals could be distributed across the network using event-based communication schemes similar to, for example, the address-event representation protocol (59).
As mentioned above, EventProp can be extended to multi-compartment neuron models as used in a recent neuromorphic architecture (60).
We used a two-layer feed-forward architecture to demonstrate learning using EventProp.The algorithm can, however, compute the gradient for arbitrary recurrent or convolutional architectures.Its computational and spatial complexity scales linearly with network size (assuming constant average firing rates per neuron), analogous to backpropagation in non-spiking artificial neural networks.The performance in more complex tasks therefore hinges on the general efficacy of gradient-based optimization in spiking networks.As mentioned above, gradients with respect to loss functions defined in terms of spike times or membrane potentials ignores the presence of critical parameters where spikes appear or disappear.We suggest that studying regularization techniques which deal with this fundamental issue in a targeted manner could enable powerful learning algorithms for spiking networks.By providing a theoretical foundation for backpropagation in spiking networks, we support future research that combines such regularization techniques with the computation of exact gradients.

Partial Derivatives in a Hybrid System
In the following, we use the example of a bouncing ball (Figure 1 A) to illustrate the calculation of partial derivatives in a dynamical system with state discontinuities.A general treatment of the topic is given in (8) or (61).The discontinuities occurring in the leaky integrate-and-fire neuron are treated analogously in our derivation of the gradient (section 4.3).
The differential equation describing the bouncing ball with height y is with gravitational acceleration g.Defining the ball's velocity as v ≡ ẏ, this is equivalent to a two-dimensional system The initial conditions are where y 0 > 0 is the parameter of interest defining the ball's initial height.The given equations determine the state trajectory y(t) up to the moment of impact with the ground at y = 0. Likewise, the trajectories of the partial derivatives with respect to y 0 are given by differentiation of eqs.( 6) and ( 7) (62), with initial conditions ∂v ∂y 0 (0) = 0, (9a) The state discontinuity occurs when the ball hits the ground and we have at the time of impact t r .The ball is inelastically reflected, losing a fraction of its energy.Specifically, the system is re-initialized as where − and + denote the state before and after the transition (v ± , y ± are functions of t r and y 0 ).Equation (10) and eq.( 11) together uniquely determine the partial derivatives after the reflection.The implicit function theorem (63) applied to eq. ( 10) guarantees (because v = 0) the existence of a function t r (y 0 ) that locally describes how the impact time changes with y 0 , with its derivative given by Likewise, the implicit function theorem applies to eq. ( 11) (because v = 0, v = 0), yielding after differentiation The partial derivatives after the transition can now be found by solving the system of equations given by eqs.(11) to (13), where we have used ÿ = −g.Equation ( 14) provides the initial conditions for the integration of the partial derivatives after the transition; subsequent ground impacts can be treated equivalently.Figure 1 A illustrates the behaviour of y(t) and ∂y ∂y0 (t) using trajectories calculated numerically using the equations given here.

Adjoint Method
We apply the adjoint method to a continuous, first order system of ordinary differential equations and refer the reader to (64; 65) for a more general setting.Consider an N -dimensional dynamical system x : t → x(t) ∈ R N with parameters p ∈ R P defined by the system of implicit first order ordinary differential equations and constant initial conditions G(x(0)) = 0 where F , G are smooth vector-valued functions.
We are interested in computing the gradient of a loss that is the integral of a smooth function l over the trajectory of x, We have where • is the dot product and the dynamics of the partial derivatives ∂x ∂pi are given by applying Gronwall's theorem (62), Computing x(t) along with ∂x ∂pi (t) using eqs.( 15) and ( 18) allows us to calculate the gradient in eq. ( 17) in a single forward pass.However, this procedure can incur prohibitive computational cost.When considering a recurrent neural network with N neurons and P = N 2 synaptic weights, computing ∂x ∂pi (t) for all parameters requires storing and integrating P N = N 3 partial derivatives.
The adjoint method allows us to avoid computing P N partial derivatives in the forward pass by instead computing N adjoint variables λ(t) in an additional backward pass.We add a Lagrange multiplier λ : t → λ(t) ∈ R N that constrains the system dynamics as given in eq. ( 15), Along trajectories where eq. ( 15) holds, λ can be chosen arbitrarily without changing L or its derivative.We get Using partial integration, we have By setting λ(T ) = 0, the boundary term vanishes because we chose parameter independent initial conditions ( ∂x ∂pi (0) = 0).The gradient becomes By choosing λ to fulfill the adjoint differential equation we are left with The gradient can therefore be computed using eq.( 24), where the adjoint state variable λ is computed from t = T to t = 0 as the solution of the adjoint differential equation eq. ( 23) with initial condition λ(T ) = 0.This corresponds to backpropagation through time (BPTT) in discrete time artificial neural networks.

Derivation of Gradient
We apply the adjoint method (section 4.2) to the case of a spiking neural network (i.e., a hybrid, discontinuous system with parameter dependent state transitions).The following derivation is specific to the model given in table 1.A fully general treatment of (adjoint) sensitivity analysis in hybrid systems can be found in ( 8) or (10).
The differential equations defining the free dynamics in implicit form are where f V , f I are again vectors of size N .We now split up the loss integral in eq. ( 1) at the spike times t post and use vectors of Lagrange multipliers λ V , λ I that fix the system dynamics f V , f I between transitions.
where we set t post 0 = 0 and t post Npost+1 = T and x • y is the dot product of two vectors x, y.Note that because f V , f I vanish along all considered trajectories, λ V and λ I can be chosen arbitrarily without changing L or its derivative.Using eq. ( 25) we have, as per Gronwall's theorem (62), where we have used the fact that the derivatives commute, ∂ where l ± V,k is the voltage-dependent loss evaluated before (−) or after (+) the transition and we have used that f V = f I = 0 along all considered trajectories.Using partial integration, we have Collecting terms in ∂V ∂wji , ∂I ∂wji , we have Since the Lagrange multipliers λ V (t), λ I (t) can be chosen arbitrarily, this form allows us to set the dynamics of the adjoint variables between transitions.Since the integration of the adjoint variables is done from t = T to t = 0 in practice (i.e., reverse in time), it is practical to transform the time derivative as d dt → − d dt .Denoting the new time derivative by , we have The integrand in eq. ( 31) therefore vanishes along the trajectory and we are left with a sum over the transitions.Since the initial conditions of V and I are assumed to be parameter independent, we have ∂V ∂wji = ∂I ∂wji = 0 at t = 0. We set the initial condition for the adjoint variables to be λ V (T ) = λ I (T ) = 0 to eliminate the boundary term for t = T .We are therefore left with a sum over transitions ξ k evaluated at the transition times t post k , with the definition We proceed by deriving the relationship between the adjoint variables before and after each transition.Since the computation of the adjoint variables happens in reverse time in practice, we provide λ − in terms of λ + .
Consider a spike caused by the nth neuron, with all other neurons m = n remaining silent.We start by first deriving the relationships between ∂V + ∂wji , ∂V − ∂wji and ∂I + ∂wji , ∂I − ∂wji .Membrane Potential Transition By considering the relations between V + , V − and V + , V − , we can derive the relation between ∂V + ∂wji and ∂V − ∂wji at each spike.Each spike at t post is triggered by a neuron's membrane potential crossing the threshold.We therefore have, at t post , (V − ) n − ϑ = 0. ( This relation defines t post as a differentiable function of w ji via the implicit function theorem (illustrated in fig.5, see also 66), under the condition that ( V − ) n = 0. Differentiation of this relation yields Since we only allow transitions for ( V − ) n = 0, we have Note that corresponding relations were previously used to derive gradient-based learning rules for spiking neuron models (67; 21; 22; 23; 27); in contrast to the suggestion in (21), eq. ( 37) is not an approximation but rather an exact relation at all non-critical parameters and invalid at all critical parameters.
Because the spiking neuron's membrane potential is reset to zero, we have This implies by differentiation Using eq. ( 37), this allows us to relate the partial derivative after the spike to the partial derivative before the spike, Since we have (V + ) m = (V − ) m for all other, non-spiking neurons m = n, it holds that Because the spiking neuron n causes the synaptic current of all neurons m = n to jump by w mn , we have and therefore get with eq. ( 36) Synaptic Current Transition The spiking neuron n causes the synaptic current of all neurons m = n to jump by the corresponding weight w mn .We therefore have By differentiation, this relation implies the consistency equations for the partial derivatives ∂I ∂wji with respect to the considered weight w ji , where δ ji is the Kronecker delta.Because we get with eq. ( 36) With (I + ) n = (I − ) n and ( İ+ ) n = ( İ− ) n , we have Using the relations of the partial derivatives from eqs. ( 37), ( 40), ( 44), ( 49) and ( 50) in the transition equation eq. ( 34), we now derive relations between the adjoint variables.Collecting terms in the partial derivatives and writing the index of the spiking neuron for the kth spike as n(k), we have This form dictates the jumps of the adjoint variables for the spiking neuron n and all other, silent neurons m,

Figure 1 :
Figure1: We derive the precise analogue to backpropagation for spiking neural networks by applying the adjoint method together with the jump conditions for partial derivatives at state discontinuities, yielding exact gradients with respect to loss functions based on membrane potentials or spike times.A, B: Dynamical systems with parameter-dependent discontinuous state transitions typically have discontinuous partial derivatives of state variables with respect to system parameters (1), as is the case for the two examples shown here.Both examples model dynamics occurring on short timescales, namely inelastic reflection and the neuronal spike mechanism, using an instantaneous state transition.We denote quantities evaluated before and after a given transition by − and +.In A, a bouncing ball starts at height y 0 > 0 and is described by ÿ = −g with gravitational acceleration g.It is inelastically reflected as ẏ+ = −0.8ẏ− as soon as y − = 0 holds, causing the partial derivative with respect to y 0 to jump as ∂y + ∂y0 = −0.8∂y − ∂y0 (see section 4.1).In B, a leaky integrate-and-fire neuron described by the system given in table 1 with initial conditions I(0) = w, V (0) = 0 resets its membrane potential as V + = 0 when V − = ϑ holds, causing the partial derivative to jump as

Figure 2 :
Figure2: Illustration of EventProp-based gradient calculation in two leaky integrate-and-fire neurons connected with weight w and a spike-time dependent loss L. The forward pass (B, C) computes the spike times for both neurons and the backward pass (D-G) backpropagates errors at spike times, yielding the gradient as given in eq.(2).A: The upper neuron receives 100 independent Poisson spike trains with frequency 200 Hz across randomly initialized weights and is connected to the lower neuron via a single weight w.The loss L is a sum of the spike times of the lower neuron.B, C: Membrane potential of upper and lower neuron.Spike times of the upper neuron are indicated using arrows.D, E: Adjoint variable λ I of upper and lower neuron.The lower neuron backpropagates its error signal λ V − λ I at the upper neuron's spike times (indicated by arrows).F, G: Accumulated gradient for one of the 100 input weights of the upper neuron and the weight w connecting the upper and lower neuron.EventProp computes the adjoint variables from t = T to t = 0 and accumulates the gradients by sampling −τ syn λ I when spikes are transmitted across the respective weight.The gradients computed in this way match the gradients computed via central differences (dashed lines) up to a relative deviation of less than 10 −7 .

Figure 3 :
Figure 3: We used EventProp and a time-to-first-spike loss function to train a two-layer leaky integrate-and-fire network on the Yin-Yang dataset.A: Illustration of the two-dimensional training dataset.The three different classes are shown in red, green and blue.This dataset was encoded using spike time latencies (see D). B, C: Training results in terms of test error and loss averaged over 10 different random seeds (individual traces shown as grey lines).D: Data points (x, y) were transformed into (x, 1 − x, y, 1 − y) and encoded using spike time latencies.We added a fixed spike at time t bias .E: Spike time latencies ∆t of the three output neurons (encoding the blue, red or green class) after training, for all samples in the test set and a specific random seed.Latencies are relative to the first spike among the three neurons and given in units of t max .A latency of zero (bright yellow dots) implies that the corresponding neuron fired the first spike, determining the class assignment.Missing spikes are denoted using green crosses.

Figure 4 :
Figure 4: We used EventProp and a two-layer network composed of a hidden layer of leaky integrate-and-fire neurons and a readout layer of non-firing neurons to classify the MNIST dataset, with the readout neuron with the largest voltage deflection determining the class assignment.A, B: Training results in terms of test error and loss averaged over 10 different random seeds (individual traces shown as grey lines).C: Confusion matrix after training for a specific random seed and using the test set.D: Voltage traces of all readout layer neurons for three different samples from the test set, where voltage traces of neurons corresponding to wrong labels are plotted using dashed lines.
weights are fixed and have no time dependence).The gradient then becomes, by application of the Leibniz integral rule,

Figure 5 :
Figure 5: In this sketch, the relation v(t, w) − ϑ = 0 defines an implicit function (black line along which dv = 0).The critical point where the gradient diverges is shown in red.

Table 3 :
Comparison of previously published classification results on the MNIST dataset for spiking neural networks that are trained using supervised learning with a single, fully connected (non-convolutional) hidden layer and temporal encoding of input data.The second column provides the number of hidden neurons.