Exact Gradients for Stochastic Spiking Neural Networks Driven by Rough Signals

Authors: christian  holberg
cristopher  salvi
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 18:28:34

Abstract: We introduce a mathematically rigorous framework based on rough path theory to model stochastic spiking neural networks (SSNNs) as stochastic differential equations with event discontinuities (Event SDEs) and driven by càdlàg rough paths. Our formalism is general enough to allow for potential jumps to be present both in the solution trajectories as well as in the driving noise. We then identify a set of sufficient conditions ensuring the existence of pathwise gradients of solution trajectories and event times with respect to the network's parameters and show how these gradients satisfy a recursive relation. Furthermore, we introduce a general-purpose loss function defined by means of a new class of signature kernels indexed on càdlàg rough paths and use it to train SSNNs as generative models. We provide an end-to-end autodifferentiable solver for Event SDEs  and make its implementation available as part of the $\texttt{diffrax}$ library. Our framework is, to our knowledge, the first enabling gradient-based training of SSNNs with noise affecting both the spike timing and the network's dynamics.

Main paper: 
 introduction: 
 Stochastic differential equations exhibiting event discontinuities (Event SDEs) and driven by noise processes with jumps are an important modelling tool in many areas of science. One of the most notable examples of such systems is that of stochastic spiking neural networks (SSNNs). Several models for neuronal dynamics have been proposed in the computational neuroscience literature with the stochastic leaky integrate-and-fire (SLIF) model being among the most popular choices [19, 56] . In its simplest form, given some continuous input current i t on [0, T ], the dynamics of a single SLIF neuron consist of an Ornstein-Uhlenbeck process describing the membrane potential as well as a threshold for spike triggering and a resetting mechanism [33] . In particular, between spikes, the dynamics of the membrane potential v t is given by the following SDEEQUATIONwhere µ > 0 is a parameter and B t is a standard Brownian motion. The neuron spikes whenever the membrane potential v hits the threshold ψ > 0 upon which v is reset to 0. Alternatively, one can model the spike times as a Poisson process with intensity λ : R → R + depending on the membrane potential v t . A common choice is λ(v) = exp((vψ)/β) [50, 26, 27, 24] .A notorious issue for calibrating Event SDEs such as SSNNs is that the implicitly defined event discontinuities, e.g., the spikes, make it difficult to define derivatives of the solution trajectories and of the event times with respect to the network's parameters using classical calculus rules. This issue is exacerbated when the dynamics are stochastic in which case the usual argument relying on the implicit function theorem, used for instance in [6, 25] , is no longer valid.In this paper, we introduce a mathematically rigorous framework to model SSNNs as SDEs with event discontinuities and driven by càdlàg rough paths, without any prior knowledge of the timing of events. The mathematical formalism we adopt is that of rough path theory [38] , a modern branch of stochastic analysis providing a robust solution theory for stochastic dynamical systems driven by noisy, possibly discontinuous, rough signals. Although Brownian motion is a prototypical example, these signals can be far more irregular (or rougher) than semimartingales [17, 16, 39] .Equipped with this formalism, we proceed to identify sufficient conditions under which the solution trajectories and the event times are differentiable with respect to the network's parameters and obtain a recursive relation for the exact pathwise gradients in Theorem 3.2. This is a strict generalization of the results presented in [6] and [25] which only deal with ordinary differential equations (ODEs). Furthermore, we define Marcus signature kernels as extensions of continuous signature kernels [52] to càdlàg rough paths and show their characteristicness. We then make use of this class of kernels indexed on discontinuous trajectories to define a general-purpose loss function enabling the training of SSNNs as generative models. We provide an end-to-end autodifferentiable solver for Event SDEs (Algorithm 1) and make its implementation available as part of the diffrax library [28] .Our framework is, to our knowledge, the first allowing for gradient-based training of a large class of SSNNs where a noise process can be present in both the spike timing and the network's dynamics. In addition, we believe this work is the first enabling the computation of exact gradients for classical SNNs whose solutions are approximated via a numerical solver (not necessarily based on a Euler scheme). In fact, previous solutions are based either on surrogate gradients [46] or follow an optimisethen-discretise approach deriving adjoint equations [56] , the latter yielding exact gradients only in the scenario where solutions are available in closed form and not approximated via a numerical solver. 1Finally, we discuss how our results lead to bioplausible learning algorithms akin to e-prop [2] . 
related work: 
 Neural stochastic differential equations (NSDEs) The intersection between differential equations and deep learning has become a topic of great interest in recent years. A neural ordinary differential equation (NODE) is an ODE of the form dy t = f θ (y t )dt started at y 0 ∈ R e using a parametric Lipschitz vector field f θ : R e → R e , usually given by a neural network [5] . Similarly, a neural stochastic differential equation (NSDE) is an SDE of the form dy t = µ θ (y t )dt + σ θ (y t )dB t driven by a d-dimensional Brownian motion B, started at y 0 ∈ R e , and with parametric vector field µ θ : R e → R e and σ θ : R e → R e×d that are Lip 1 and Lip 2+ϵ continuous respectively 2 . Rough path theory offers a way of treating ODEs, SDEs, and more generally differential equations driven by signals or arbitrary (ir)regularity, under the unified framework of rough differential equations (RDEs) [44, 22] . For an account on applications of rough path theory to machine learning see [4, 14, 51] .Training techniques for NSDEs Training a NSDE amounts to minimising over model parameters an appropriate notion of statistical divergence between a distribution of continuous trajectories generated by the NSDE and an empirical distribution of observed sample paths. Several approaches have been proposed in the literature, differing mostly in the choice of discriminating divergence. SDE-GANs, introduced in [29] , use the 1-Wasserstein distance to train a NSDE as a Wasserstein-GAN [1] .Latent SDEs [36] train a NSDE with respect to the KL divergence via variational inference and can be interpreted as variational autoencoders. In [23] the authors propose to train NSDEs non-adversarially using a class of maximum mean discrepancies (MMD) endowed with signature kernels [30, 52] . Signature kernels are a class of characteristic kernels indexed on continuous paths that have received increased attention in recent years thanks to their efficiency for handling path-dependent problems [35, 54, 10, 53, 9, 48, 41] . For a treatment of this topic we refer the interested reader to [4, Chapter 2] . These kernels are not applicable to sample trajectories of SSNNs because of the lack of continuity.Backpropagation through NSDEs Once a choice of discriminator has been made, training NSDEs amounts to perform backpropagation through the SDE solver. There are several ways to do this. The first option is simply to backpropagate through the solver's internal operations. This method is known as discretise-then-optimise; it is generally speaking fast to evaluate and produces accurate gradients, but it is memory-inefficient, as every internal operation of the solver must be recorded. A second approach, known as optimise-then-discretise, computes gradients by deriving a backwards-in-time differential equation, the adjoint equation, which is then solved numerically by another call to the solver. Not storing intermediate quantities during the forward pass enables model training at a memory cost that is constant in depth. Nonetheless, this approach produces less accurate gradients and is usually slower to evaluate because it requires recalculating the forward solutions to perform the backward pass. A third way of backpropagating through NDEs is given by algebraically reversible solvers, offering both memory and accuracy efficiency. We refer to [28] for further details.Differential equations with events Many systems are not adequately modelled by continuous differential equations because they experience jump discontinuities triggered by the internal state of the system. Examples include a bouncing ball or spiking neurons. Such systems are often referred to as (stochastic) hybrid systems [21, 37] . When the differential equation is an ODE, there is a rich literature on sensitivity analysis aimed at computing derivatives using the implicit function theorem [11, 12] . If, additionally, the vector fields describing the hybrid system are neural networks, [6] show that NODEs solved up until first event time can be implemented as autodifferentiable blocks and [25] derive the corresponding adjoint equations. Nonetheless, none of these works cover the more general setting of SDEs. The only work, we are familiar with, dealing with sensitivity analysis in this setting is [47] , although focused on the problem of optimal control.Training techniques for SNNs Roughly speaking, these works can be divided into two strands. The first, usually referred to as backpropagation through time (BPTT), starts with a Euler approximation of the SNN and does backpropagation by unrolling the computational graph over time; it then uses surrogate gradients as smooth approximations of the gradients of the non-differentiable terms. [58, 46, 40] . This approach is essentially analogous to discretise-then-optimise where the backward pass uses custom gradients for the non-differentiable terms. The second strand computes exact gradients of the spike times using the implicit function theorem. These results are equivalent to optimise-then-discretise and can be used to define adjoint equations as in [56] or to derive forward sensitivities [34] . However, we note that, unless solution trajectories and spike times of the SNN are computed exactly, neither method provides the actual gradients of the implemented solver. Furthermore, the BPTT surrogate gradient approach only covers the Euler approximation whereas many auto-differentiable differential equation solvers are available nowadays, e.g. in diffrax.Finally, there is a lot of interest in developing bioplausible learning algorithms where weights can be updated locally and in an online fashion. Notable advances in this direction include [2, 57] . To the best of our knowledge, none of these works cover the case of stochastic SNNs where the neuronal dynamics are modeled as SDEs instead of ODEs.Question: Does the paper discuss the limitations of the work performed by the authors? 
stochastic spiking neural networks as event sdes: 
 We shall in this paper be concerned with SDEs where solution trajectories experience jumps triggered by implicitly defined events, dubbed Event SDEs. The prototypical example that we come back to throughout is the SNN model composed of SLIF neurons. Here the randomness appears both in the inter-spike dynamics as well as in the firing mechanism. To motivate the general definitions and concepts we start with an informal introduction of SSNNs.To achieve a more bioplausible model of neuronal behaviour, one can extend the simple deterministic LIF model by adding two types of noise: a diffusion term in the differential equation describing inter-spike behaviour [33] and stochastic firing [50, 27] . That is, the potential is modelled by eq. ( 1). Instead of firing exactly when the membrane potential hits a set threshold, we model the spike times (event times) by an inhomogenous Poisson process with intensity λ : R e → R + which is assumed to be bounded by some constant C > 0. This can be phrased as an Event SDE (note that this is essentially the reparameterisation trick) by introducing the additional state variable s t satisfyingds t = λ(v t-)dt, s 0 = log uwhere u ∼ Unif(0, 1). The neuron spikes whenever s t hits 0 from below at which point the membrane potential is reset to a resting level and we sample a new initial condition for s t . We can denote this first spike time by τ 1 and repeat the procedure to generate a sequence of spike times τ 1 < τ 2 < ... In practice, we reinitialize s t at log uα for some α > 0. It can then be shown thatP (t < τ n+1 |F τn ) = min 1, exp α - t τn λ(v t-)dt for t ∈ [τ n , τ n+1 ).It follows that τ n+1τ n ≥ α/C a.s., i.e. α controls the refractory period after spikes, a large value indicating a long resting period.We can then build a SSNN by connecting such SLIF neurons in a network. In particular, apart from the membrane potential, we now also model the input current of each neuron as affected by the other neurons in the network. Let K ≥ 1 denote the total number of neurons. We model neuron k ∈ [K] be the three dimensional vector y k = (v k , i k , s k ) the dynamic of which in between spikes is given byEQUATIONwhereB k is a standard two-dimensional Brownian motion, σ = (σ 1 , σ 2 ) ∈ R 2×2 , µ = (µ 1 , µ 2 ) ∈ R 2, and λ(•; ξ) : R → R + is an intensity function. As before, neuron k fires (or spikes) whenever s k hits zero from below. Apart from resetting the membrane potential, this event also causes spikes to propagate through the network in a such a way that a spike in neuron k will increment the input current of neuron j by w kj . Here w ∈ R K×K is a matrix of weights representing the synaptic weights in the neural network. If one is only interested in specific network architectures such as, e.g., feed-forward, this can be achieved by fixing the appropriate entries in w at 0.As presented here, there is no way to model information coming into the network. But this would only require a minor change. Indeed, by adding a suitable control term to eq. ( 2) we can model all relevant scenarios. Since this does not change the theory in any meaningful way (the general theory in Appendix B covers RDEs so an extra smooth control is no issue), we only discuss the more simple model given without any additional input currents.Definition 3.1 (Event SDE). Let N ∈ N be the number of events. Let y 0 ∈ R e be an initial condition. Let µ : R e → R e and σ : R e → R e×d be the drift and diffusion vector fields. Let E : R e → R and T : R e → R e be an event and transition function respectively. We say that y, (τ n ) N n=1 is a solution to the Event SDE parameterised by (y 0 , µ, σ, E, T , N ) if y T = y N T ,EQUATIONwith E(y n τn ) ̸ = 0 andEQUATIONEQUATIONwhere B t is a d-dimensional Brownian motion and (4), ( 5) are Stratonovich SDEs.In words, we initialize the system at y 0 , evolve it using (4) until the first time τ 1 at which an event happens E(y 0 τ1 ) = 0. We then transition the system according to y 1 τ1 = T y 0 τ1-and evolve it according to (5) until the next event is triggered. We note that Definition 3.1 can be generalised to multiple event and transition functions. Also, the transition function can be randomised by allowing it to have an extra argument u ∼ Unif([0, 1]). As part of the definition we require that there are only finitely many events and that an event is not immediately triggered upon transitioning.Existence of strong solutions to Event SDEs driven by continuous semimartingales has been studied in [31, Theorem 5.2 ] and [32] . Under sufficient regularity of µ and σ, a unique solution to (4) exists. We need the following additional assumptions: Assumption 3.1. There exists c > 0 such that for all s ∈ (0, T ) and a ∈ im T it holds that inf{t > s : E(y t ) = 0} > c where y t is the solution to 4 started at y s = a Assumption 3.2. It holds that T (ker E) ∩ E = ∅. Assumptions 3.1 and 3.2 simply ensure that an event cannot be triggered immediately upon transitioning. This holds in most settings of interest. For example, for the usual deterministic LIF neuron im T = 0 and ker E = 1 and the duration of the refractory period is directly linked to c in Assumption 3.1. Theorem 3.1 (Theorem 5.2, [31] ). Under Assumptions 3.1-3.2 and with µ ∈ Lip 1 and σ ∈ Lip γ for γ > 2, there exists a unique solution (y, (τ n ) N n=1 ) to the Event SDE of Definition 3.1.The definitions and results of this section can be extended to differential equations driven by random rough paths, and in particular, to cases where the driving noise exhibits jumps. In the latter case, it is important to note that the resulting Event SDE will exhibit two types of jumps: the ones given apriori by the driving noise and the ones that are implicitly defined through the solution (what we call events). In fact, we develop the main theory of Event RDEs in Appendix A in the more general setting of RDEs driven by càdlàg rough paths. The rough path formalism enables a unified treatment of differential equations driven by noise signals of arbitrary (ir-)regularity, and makes all proofs simple and systematic. In particular, it allows us to handle cases where the diffusion term is driven by a finite activity Lévy process (e.g, a homogeneous Poisson processes highly relevant in the context of SNNs).We are interested in optimizing a continuously differentiable loss function L whose input is the solution of a parameterised Event SDE. As for Neural ODEs, the vector fields, µ, σ, and the event and transition functions E, T , might depend on some learnable parameters θ. We can move the parameters θ of the Event RDE inside the initial condition y 0 by augmenting the dynamics with the additional state variable θ t satisfying dθ t = 0 and θ 0 = θ. Thus, as long as we can compute gradients with respect to y 0 , these will include gradients with respect to such parameters. We then require the gradients ∂ y0 L, if they exist. For this, we need to be able to compute the Jacobians ∂y n t := ∂ y0 y n t of the inter-event flows associated to the dynamics of y n t and the derivatives ∂τ n := ∂ y0 τ n . We assume that the event and transition functions E and T are continuously differentiable. Apriori, it is not clear under what conditions such quantities exist and even less how to compute them. This shall be the focus of the present section. We will need the following running assumptions. Assumption 3.3. σ(T (y)) -∇T (y)σ(y) = 0 for all y ∈ ker E. Assumption 3.4. ∇E(y)σ(y) = 0 for all y ∈ ker E. Assumption 3.5. ∇E(y)µ(y) ̸ = 0 for all y ∈ ker E. Assumption 3.4 and 3.5 ensure that the event times are differentiable. Intuitively, they state that the event condition is hit only by the drift part of the solution. Assumption 3.4 holds for example if the event functions depend only on a smooth part of the system. Assumption 3.3 is what allows us to differentiate through the event transitions. Theorem 3.2. Let Assumptions 3.1-3.5 be satisfied and (y, (τ n ) N n=1 ) the solution to the Event SDE parameterized by (y 0 , µ, σ, E, T , N ). Then, almost surely, for any n ∈ [N ], the derivatives ∂τ n and the Jacobians ∂y n t exist and admit the following recursive expressionsEQUATIONEQUATION)where ∂y n t and ∂τ n are the total derivatives of y t n and τ n with respect to the initial condition y 0 , ∂ y n τn y n t denotes the partial derivative of the flow map of eq. (5) with respect to its initial condition, and ∇T ∈ R e×e and ∇E ∈ R 1×e are the Jacobians of T and E. Remark 3.1. If the diffusion term is absent we recover the gradients in [6] . In this case, the assumptions of the theorem are trivially satisfied. Note however, that the result, as stated here, is slightly different since we are considering repeated events. Remark 3.2. The recursive nature of ( 6) - (7) suggest a way to update gradients in an online fashion by computing the forward sensitivity along with the state of the Event SDE. In traditional machine learning applications (e.g. NDEs) forward mode automatic differentiation is usually avoided due to the fact that the output dimension tends to be orders of magnitude smaller than the number of parameters [28] . However, for (S)SNNs this issue can be partly avoided as discussed in Section 4.4.Returning now to the SSNN model introduced in Section 3.1 we find that it is an Event SDE with K different event functions given by E k (y) = s k and corresponding transition functions given byT k (y) = T 1 k (y 1 ), . . . , T K k (y K )whereT j k (y j ) = (v j , i j + w kj , s j ) if j ̸ = k and T k k (y k ) = (v k -v reset , i k , log u -α)where v reset > 0 is a constant determining by what amount the membrane potential is reset. The addition of the constant α > 0 controlling the refractory period ensures that Assumption 3.2 and 3.2 are satisfied. Stochastic firing smooths out the event triggering so that Assumption 3.5 and 3.4 hold. Finally, one can check that the combination of constant diffusion terms and the given transition functions satisfies Assumptions 3.3. Note that setting v k t exactly to 0 upon spiking would break Assumption 3.3. If one is interested in such a resetting mechanism it suffices to pick a diffusion term σ 1 (y k ) that satisfies σ(0) = 0. To sum up, solutions (in the sense of Def. 3.1) of the SSNNs exist and are unique. In addition, the trajectories and spike times are almost surely differentiable satisfying ( 6) and (7).Theorem 3.2 gives an expression for the gradients of the event times as well as the Event SDE solution. In practice, analytical expressions for gradients are often not available and one has to resort to numerical solvers. Three solutions suggest themselves:1. There are multiple autodifferentiable differential equation solvers (such as diffrax [28] ) that provide differentiable numerical approximations of the flows ∂ y n τn y n t . We shall write SDESolve(y 0 , µ, σ, s, t) for a generic choice of such a solver. Furthermore, if RootFind(y 0 , f ) is a differentiable root finding algorithm (here f : (y, t) → R should be differentiable in both arguments and RootFind(y 0 , f ) returns t * ∈ R such that f (y 0 , t * ) = 0), then we can define a differentiable map E : y 0 → y * by t * = RootFind (y 0 , E(SDESolve (•, µ, σ, s, •))) , y * = SDESolve (y 0 , µ, σ, s, t * ) .Consequently, EventSDESolve(y 0 , µ, σ, E, T , N ) can be implemented as subsequent compositions of T • E (see Algorithm 1). This is a discretise-then-optimise approach [28] .2. Alternatively, one can use the formulas (6) and (7) directly as a replacement of the derivatives. This is the approach taken in e.g. [6] . To be precise, one would replace all the derivatives of the flow map (terms of the sort ∂ y n τn y n t ) with the derivatives of the given numerical solver. This approach is a solution between discretise-then-optimise and optimise-then-discretise. 3. Finally, one could apply the adjoint method (or optimise-then-discretise) as done for deterministic SNNs in [56] by deriving the adjoint equations. These adjoint equations define another SDE with jumps which is solved backwards in time. Between events the dynamics are exactly as in the continuous case so one just needs to specify the jumps of the adjoint process. This can be done by referring to ( 6) and (7) .Remark 3.3. One thing to be careful of with the discretise-then-optimise approach is that the SDE solver will compute time derivatives in the backward pass, although the modelled process is not time differentiable. Assumptions 3.4 and 3.3 should in principle guarantee that these derivatives cancel out (see Appendix B), yet this might not necessarily happen at the level of the numerical solver because of precision issues. This is essentially due to the fact that approximate solutions provided by numerical solvers are in general not flows. Thus, when the path driving the diffusion term is very irregular, the gradients can become unstable. In practice we found this could be fixed by setting the gradient with respect to time of the driving Brownian motion to 0 and picking a step size sufficiently small. Remark 3.4. In the context of SNNs, Algorithm 1 is actually a version of exact backpropagation through time (BPTT) of the unrolled numerical solution. Contrary to popular belief, this illustrates that one can compute exact gradients of numerical approximations of SNNs without the need to resort to surrogate gradient functions. Of course, this does not alleviate the so-called dead neuron problem. However, this ceases to be a problem when stochastic firing is introduced. In fact, surrogate gradients can be related to stochastic firing mechanisms and expected gradients [20] .Remark 3.5. One the one hand, the EventSDESolve algorithm as presented here scales poorly in the number of events since it requires doing a full SDESolve and an additional RootFind each time an event occurs. This problem becomes especially prevalent for SSNNs with a large number of neurons since in this case an event is triggered every time a single neuron spikes and the inter-spike SDE that needs to be solved is high-dimensional. On the other hand, there are multiple ways to mitigate this issue. Firstly, one could relax the root-finding step and simply trigger a spike as soon as e ≥ 0 and take this as the spike time. For the backward pass one could then solve the adjoint equations (for which you need need to store the spike times in the forward pass). The resulting algorithm would be similar to the one derived in [55] for deterministic SNNs. Secondly, for special architectures such as a feed-forward network, given the spikes from the previous layer, one could solve the EventSDE for each neuron in the current layer independently of all other neurons. This would imply that a forward (or backward) pass of the entire SSNN scales as O(KS) where S is the cost of the forward (or backward) pass of a single neuron and K is the number of neurons.Input y 0 , µ, σ, E, T , N, t 0 , ∆t, T 1: y ← y 0 2: n ← 0 3: e ← E(y) 4: while n < N and t 0 < T do 5:while e < 0 do ▷ We assume for simplicity that e ≤ 0 6:y 0 ← y 7:y ← SDESolveStep(y 0 , µ, σ, t 0 , ∆t)8: t 0 ← t 0 + ∆t 9:e ← E(y) ▷ Update value of event function 10:end while To train SSNNs we will adopt a similar technique as in [23] , where the authors propose to train NSDEs non-adversarially using a class of maximum mean discrepancies (MMD) endowed with signature kernels [52] indexed on spaces of continuous paths as discriminators. However, as we mentioned in the introduction, classical signature kernels are not directly applicable to the setting of SSNNs as the solution trajectories not continuous. To remedy this issue, in Appendix C, we generalise signature kernels to Marcus signature kernels indexed on discontinuous (or càdlàg) paths. We note that our numerical experiments only concern learning from spike trains, which are càdlàg paths of bounded variation. Yet, the Marcus signature kernel defined in Appendix C can handle more general càdlàg rough paths.11: t * n+1 ← RootFind (y 0 , E(SDESolveStep (•, µ, σ, t 0 -∆t, •))) ▷ Find exact event time 12: y * n+1 ← SDESolveStep(y 0 , µ, σ, t 0 -∆ t , t * n+1 ) ▷ ComputeThe main idea goes as follows. If x is a càdlàg path, one can define the Marcus signature S(x) in the spirit of Marcus SDEs [42, 43] as the signature of the Marcus interpolation of x. The general construction is given in Appendix A. The Marcus signature kernel is defined as the inner product k(x, y) = ⟨S(x), S(y)⟩ of Marcus signatures S(x), S(y) of two càdlàg paths x, y. As stated in the first part of Theorem C.1, this kernel is characteristic on regular Borel measures supported on compact sets of càdlàg paths. In particular, this implies that the resulting maximum mean discrepancy (MMD)d k (µ, ν) 2 = E x,x ′ ∼µ k(x, x ′ ) -2E x,y∼µ×ν k(x, x ′ ) + E y,y ′ ∼ν k(y, y ′ )satisfies the property d k (µ, ν) 2 = 0 ⇐⇒ µ = ν for any two compactly supported measures µ, ν.Nonetheless, characteristicness ceases to hold when one considers measures on càdlàg paths that are not compactly supported. In [8] the authors address this issue for continuous paths by using the so-called robust signature. They introduce a tensor normalization Λ ensuring that the range of the robust signature Λ • S remains bounded. The robust signature kernel is then defined as the inner product k Λ (x, y) = ⟨Λ • S(x), Λ • S(y)⟩. This normalization can be applied analogously to the Marcus signature resulting in a robust Marcus signature kernel. In the second part of Theorem C.1, we prove characteristicness of k Λ for possibly non-compactly supported Borel measures on càdlàg paths. The resulting MMD is denoted by d kΛ .There are several ways of evaluating signature kernels. The most naive is to simply truncate the signatures at some finite level and then take their inner product. Another amounts to solve a pathdependent wave equation [52] . Our experiments are compatible with both of these methods.Given a collection of observed càdlàg trajectories {x i } m i=1 ∼ µ true sampled from an underlying unknown target measure µ true , we can train an Event SDE by matching the generated càdlàg trajectories {y i } n i=1 ∼ µ θ using an unbiased empirical estimator of d k (or d kΛ ), i.e. minimising over the parameters θ of the Event SDE the following loss functionL = 1 m(m -1) j̸ =i k(x i , x j ) - 2 mn i,j k(x i , y j ) + 1 n(n -1) j̸ =i k(y i , y j ).In the context of SSNNs, the observed and generated trajectories x i 's and y i 's correspond to spike trains, which are càdlàg paths of bounded variation.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]Justification: Rigorous proofs are provided in Appendix B and the assumptions are an integral part of the paper.Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
input current estimation: 
 The first example is the simple problem of estimating the constant input current c > 0 based on a sample of spike trains in the single SLIF neuron model,dv t = µ(c -v t )dt + σdB t , ds t = λ(v t )dt,where λ(v) = exp(5(v -1), µ = 15 and σ varies. Throughout we fix the true c = 1.5 and set v reset = 1.4 and α = 0.03. We run stochastic gradient descent for 1500 steps for two choices of the diffusion constant σ. The loss function is the signature kernel MMD between a simulated batch and the sample of spike trains. 3 . The test loss is the mean absolute error between the first three average spike times. Results are given in Fig. 1 . For additional details regarding the experiments, we refer to Appendix E. In all cases backpropagation through Algorithm 1 is able to learn the underlying input current after around 600 steps up to a small estimation error. In particular, the convergence is fastest for the largest sample size and the true c is recovered for both levels of noise.Next we consider the problem of estimating the weight matrix in a feed-forward SSNN with input dimension 4, 1 hidden layer of dimension 16, and output dimension 2. The rest of the parameters are fixed throughout. We run stochastic gradient descent for 1500 steps with a batch size of 128 and for a sample size of 256, 512, and 1024 respectively. Learning rate is decreased from 0.003 to 0.001 after 1000 steps. The results are given in Fig. 2 in Appendix E. For a sample size of 512 and 1024 we are able to reach a test loss of practically 0, that is, samples from the learned model and the underlying model are more or less indistinguishable. Also, in all cases the estimated weight matrix approaches the true weight matrix. Interestingly, for the largest sample size, the model reaches the same test loss as the model trained on a sample size of 512, but their estimated weight matrices differ significantly.In the case of SSNNs, equations ( 6)-( 7) lead to a formula for the forward sensitivity where any propagation of gradients between neurons only happens at spike times and only between connected neurons (see Proposition D.1). Since the forward sensitivities are computed forward in time together with the solution of the SNN, gradients can be updated online as new data appears. As a result, between spikes of pre-synaptic neurons, we can update the gradient flow of the membrane potential and input current of each neuron using information exclusively from that neuron. For general network structures and loss functions, however, this implies that each neuron needs to store on the order of K 2 gradient flows (one for each weight in the network).On the other hand, if the adjacency matrix of the weight matrix forms a directed acyclic graph (DAG), three-factor Hebbian learning rules like those in [57, 2] are easily derived from Proposition D.1. For simplicity, consider the SNN consisting of deterministic LIF neurons and let N k t denote the spike train of neuron k, i.e., N k t is càdlàg path equal to the number of spikes of neuron k at time t. We let τ k (t) (or τ k for short) denote the last spike of neuron k before time t. We shall assume that the instantaneous loss function L t depends only on the most recent spike times τ 1 , . . . , τ K . Then,∂ w jk L t = ∂ τ k L t a jk τ k µ 1 (v k τ k -i k τ k )where a jk t is the eligibility trace and the first term can be viewed as a global modulator, that is, a top-down learning signal propagating the error from the output neurons. 4 The eligibility trace satisfiesda jk t = µ 1 b jk t -a jk t dt + v reset a jk t µ 1 (i k t -v k t ) dN k t , db jk t = -µ 2 b jk t + dN j t ,where the dN terms are to be understood in the Riemann-Stieltjes sense. In other words, the eligibility trace can be updated exclusively from the activity of the pre-and post-synaptic neurons. We note the similarity to the results derived in [2] only our result gives the exact gradients with no need to introduce surrogate gradient functions. A similar equation for deterministic SNNs was derived in [49] (see, in particular, Chapter 5). For general network structures one can use the eligibility traces as proxies for the the true derivatives ∂ wij τ k .Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: We clearly describe the main algorithm used in both experiments (Algorithm 1) as well as the experimental setups. Additionally, the code is provided in the supplementary material with clear instructions on how to setup the correct environment. 
conclusion: 
 We introduced a mathematical framework based on rough path theory to model SSNNs as SDEs exhibiting event discontinuities and driven by càdlàg rough paths. After identifying sufficient conditions for differentiability of solution trajectories and event times, we obtained a recursive relation for the pathwise gradients in Theorem 3.2, generalising the results presented in [6] and [25] which only deal with the case of ODEs. Next, we introduced Marcus signature kernels as extensions of continuous signature kernels from [52] to càdlàg rough paths and used them to define a general-purpose loss function on the space of càdlàg rough paths to train SSNNs where noise is present in both the spike timing and the network's dynamics. Based on these results, we also provided an end-to-end autodifferentiable solver for SDEs with event discontinuities (Algorithm 1) and made its implementation available as part of the diffrax repository. Finally, we discussed how our results lead to bioplausible learning algorithms akin to e-prop [2] but in the context of spike time gradients.The primary objective of the paper was to lay out the theoretical foundations of gradient-based learning with stochastic SNNs. Although we provided an initial implementation, which is well-suited for low dimensional examples, a robust version that scales to a high number of neurons is beyond the scope of the paper. Examples that require a much higher number of neurons than the two examples already discussed will be hard to handle with the discretize-then-optimize approach for the reasons given in Remark 3.5.We think there are still many interesting research directions left to explore. For instance, it would be of interest to implement the adjoint equations or to use reversible solvers and compare the results. Similarly, since our Algorithm 1 differs from the usual approach with surrogate gradients even in the deterministic setting, questions remain on how these methods compare for training SNNs. Furthermore, it would be interesting to understand to what extent the inclusion of different types of driving noises in the dynamics of SSNNs would be beneficial for learning tasks compared to deterministic SNNs. Finally, it remains to be seen whether the discussion in Section 4.4 could lead to a bio-plausible learning algorithm with comparable performance to state-of-the-art backpropagation methods and implementable on neuromorphic hardware.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: Code is provided in the supplementary material. The two main experiments can be reproduced by running the two notebooks in the snnax folder as referenced in Appendix E.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). • The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips.cc/public/guides/CodeSubmissionPolicy) for more details. • The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. • The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimiser, etc.) necessary to understand the results? Answer: [Yes] Justification: We provide the most important aspects in Section 4. The full details are given in Appendix E.• We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: We have followed the citation guidelines for referencing diffrax which is the main asset used for the current work. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [NA] Justification: Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)? Answer: [NA] Justification: Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: The appendix is structured as follows. Section A covers the basic concepts of càdlàg rough paths based on [7] extended with a few of our own definitions and results. It culminates with the definition of Event RDEs which can be viewed as generalizations of Event SDEs. Section B covers the proof of the main result, Theorem 3.2, but in the setting of Event RDEs as well as some preliminary technical lemmas needed for the proof. Section C gives a brief overview of the main concepts in kernel learning and presents our results on Marcus signature kernels along with their proofs. Section D derives the forward sensitivities of a SSNN. Finally, Section E covers all the technical details of the simulation experiments that were not discussed in the main body of the paper.Marcus integration developed in [7] preserves the chain rule and thus serves as an analog to Stratonovich integration for semi-martingales with jump discontinuities. In particular, it allows to define a canonical lift under which càdlàg semi-martingales are a.s. geometric rough paths and many of the results from the continuous case, such as universal limit theorems and stability results, carry over under suitably defined metrics.We briefly review some of the important concepts here by following the same setup as in [7] .Let C([0, T ], E) and D([0, T ], E) be the space of continuous and càdlàg paths respectively on [0, T ] with values in a metric space (E, d). For p ≥ 1, let C p ([0, T ], E) and D p ([0, T ], E) be the corresponding subspaces of paths with finite p-variation. For anyN ≥ 1, Let G N (R d ) be the step-N free nilpotent Lie group over R d endowed with the Carnot-Carathéodory metric d. Let Ω C p (R d ) := C p ([0, T ], G ⌊p⌋ (R d )) and Ω D p (R d ) := D p ([0, T ], G ⌊p⌋ (R d ))be the space of weakly geometric continuous and càdlàg p-rough paths respectively with the homogeneous p-variation metricd p (x, y) = max 1≤k≤⌊p⌋ sup D⊂[0,T ] D d(x ti,ti+1 , y ti,ti+1 ) p k k p.Define the log-linear path functionϕ : G N (R d ) × G N (R d ) → C([0, 1], G N (R d )) (a, b) → exp((1 -•) log a + • log b).where log and exp are the (truncated) tensor logarithm and exponential maps onG N (R d ). If N = 1, then G N (R d ) ∼ = R ⊕ R d and ϕ(a, b) t = (1, (1 -t)a + tb) is a straight line connecting a to b in unit time. For any x ∈ D([0, T ], G N (R d )) we can construct a continuous path x ∈ C([0, T ], G N (R d ))by adding fictitious time and interpolating through the jumps using the log-linear path function according to the following definition.Definition A.1 (Marcus interpolation). Let N ≥ 1. For x ∈ D([0, T ], G N (R d )), let τ 1 , τ 2 , . . . , τ m be the jump times of x ordered such that d(x τ1-, x τ1 ) ≥ d(x τ2-, x τ2 ) ≥ • • • ≥ d(x τm-, x τm ), where 0 ≤ m ≤ ∞ is the number of jumps. Let (r k ) be a sequence of positive scalars r k > 0 such that r = m k=1 r k < +∞. Define the discontinuous reparameterisation η : [0, T ] → [0, T + r] by η(t) = t + m k=1 r k 1 {τ k ≤t} . The Marcus augmentation x M ∈ C([0, T + r], G N (R d )) of x is the path x M s = x t , if s = η(t) for some t ∈ [0, T ], ϕ(x τ k -, x τ k ) (s-η(τ k -))/r k , if s ∈ [η(τ k -), η(τ k )) for 1 ≤ k < m + 1. The Marcus interpolation x ∈ C([0, T ], G N (R d )) of x is the path x = x M • η r where η r (t) = t(T + r)/T is a reparameterisation from [0, T ] to [0, T + r]. We can recover x from x via x = x • η x by considering the reparameterisation η x = η -1 r • η.Once the Marcus interpolation is defined we can state what we mean by a solution to a differential equation driven by a geometric càdlàg rough path. ∈ Ω C p (R d ) dŷ t = f (ŷ t )dx t , ŷ0 = a.Define the solution y ∈ D p ([0, T ], R e ) to the Marcus RDEEQUATIONto be y = ŷ • η x , where η x is the reparameterisation introduced in Definition A.1.Chevyrev and Friz [7] introduce a metric α p on Ω D p (R d ) with respect to which 1) geometric càdlàg rough paths can be approximated with a sequence of continuous paths [7, Section 3.2] and 2) the solution map (y 0 , x) → (x, y) of the Marcus RDE ( 8) is locally Lipschitz continuous [7, Theorem 3.13].We write Λ for the set of increasing bijections from [0, T ] to itself. For a λ ∈ Λ we let |λ| = sup t∈[0,T ] |λ(t) -t|. We first define the Skorokhod metric as well as a Skorokhod version of the usual p-variation metric. Definition A.3. For p ≥ 1 and x, y ∈ D p ([0, T ], E), we defineσ ∞ (x, y) = inf λ∈Λ max |λ|, sup t∈[0,T ] d((x • λ) t , y t ) , σ p (x, y) = inf λ∈Λ max {|λ|, d p (x • λ, y)} .It turns out that the topology induced by σ p is too strong. In particular, it is not possible to approximate paths with jump discontinuities with a sequence of continuous paths (see Section 3.2 in [7] ). For x ∈ Ω D p (R d ) and f = (f 1 , . . . , f d ) a family of vector fields in Lip γ-1 (R e ) with γ > p, let Φ f (y, s, t; x) denote the solution to the Marcus RDE dy t = f (y t ) ⋄ dx t initialized at y s = y and evaluated at time t. We define the setJ f = ((a, b), (a ′ , b ′ )) | a, a ′ ∈ G ⌊p⌋ (R e ), Φ f (b, 0, 1; ϕ(a, a ′ )) = b ′ .and, on it, the path functionϕ f ((a, b), (a ′ , b ′ )) t = (ϕ(a, a ′ ), Φ f (b, 0, 1; ϕ(a, a ′ ) t )) . Finally, we let D f p ([0, T ], G ⌊p⌋ (R d ) × R e) be the space of càdlàg paths z = (x, y) on G ⌊p⌋ (R d ) × R e of bounded p-variation such that (z t -, z t ) ∈ J f for all jump times t of z. To keep notation simple, we shall write D f p when this does not cause any confusion. Naturally, if y is the solution to the Marcus RDE dy t = f (y t ) ⋄ dx t , we have (x, y) ∈ D f p . For a z = (x, y) ∈ D f p we may define the Marcus interpolation by interpolating the jumps using ϕ f . Let ẑδ denote this interpolation but with r k replaced by δr k for δ > 0 and similarly for xδ with x ∈ Ω D p (R d ). Definition A.4. For f = (f 1 , . . . , f d ) a family of vector fields in Lip γ-1 (R e ) with γ > p, let z, z ′ ∈ D f p with z = (x, y) and z ′ = (x ′ , y ′ ) and defineα p (x, x ′ ) = lim δ→0 σ p (x δ , x′ δ ), α p (z, z ′ ) = lim δ→0 σ p (ẑ δ , ẑ′ δ ).Remark A.1. It is proven in [7] that in both cases the limit in α p exists, is independent of the choice of r k , and that it is indeed a metric on Ω D p (R d ) resp. D f p . Theorem A.1 (Theorem 3.13 + Proposition 3.18, [7] ). Let f = (f 1 , . . . , f d ) be a family of vector fields in Lip γ-1 (R e ) with γ > p. Then, 1. The solution mapR e × (Ω D p (R d ), α p ) → (D f p , α p ) (y 0 , x) → z = (x, y)of the Marcus RDE dy t = f (y t ) ⋄ dx t initialized at y 0 ∈ R e is locally Lipschitz.2. On sets of bounded p-variation, the solution mapR e × (Ω D p (R d ), σ ∞ ) → (D p ([0, T ], R e ), σ ∞ ) (y 0 , x) → y of the Marcus RDE dy t = f (y t ) ⋄ dx t initialized at y 0 ∈ R e is continuous. Now, let C 1 0 (R d) be the space of absolutely continuous functions on R d . Definition A.5. We define the space of geometric càdlàg p-rough paths Ω D 0,p (R d ) as the closure ofC 1 0 (R d ) in Ω D p (R d ) under the metric α p . Remark A.2. A càdlàg semi-martingale x ∈ D p ([0, T ], R d) can be canonically lifted to a geometric càdlàg p-rough path, with p ∈ [2, 3), by enhancing it with its two-fold iterated Marcus integrals, i.e.x s,t = (1, x s,t , s,t (x s -x u ) ⊗ ⋄dx u ) ∈ G 2 (R d )where the integral is defined in a similar spirit to Definition A.2 (see, for example, [18] for more information). The solution to the corresponding Marcus RDE agrees a.s. with the solution to the usual càdlàg Marcus SDE which, in turn, if x has a.s. continuous sample paths, agrees a.s. with the solution to the Stratonovich SDE. See, e.g., Proposition 4.16 in [7] .The extended tensor algebra over R d is given byT R d = ∞ n=0 R d ⊗nequipped with the usual addition + and tensor multiplication ⊗. An element a ∈ T R d is a formal series of tensors a = (a 0 , a 1 , . . . ) such that a n ∈ (R d ) ⊗n . We define the projections π n : T R d → (R d ) ⊗n given by π n (a) = a n . Let T ((R d )) be the subset of T ((R d )) such that the π 0 (a) = 1 for all a ∈ T ((R d )). Finally, we define the set of group-like elements,EQUATION) Remark A.3. Uniqueness and existence of the signature follow from the continuous analog. Indeed, by definition, ( 9) is equivalent to a continuous linear RDE. Remark A.4. The signature, as defined here, is also known as the minimal jump extension of x and was first introduced in [15] . It was further explored in [13] where it was also shown that it acts as a universal feature map.In the continuous case, it is well known that the signature characterizes paths up to tree-like equivalence. Two continuous paths x, y are said to be tree-like equivalent if there exists a continuous non-negative map h :[0, T ] → R + such that h(0) = h(T ) and ∥x s,t -y s,t ∥ ≤ h(s) + h(t) -2 inf u∈[s,t] h(u).This can be generalized to càdlàg paths in the following way. We say that two càdlàg paths x, y are tree-like equivalent, or x ∼ t y, if their their Marcus interpolations (see Def. A.1), x and ŷ, are tree-like equivalent. It is straightforward to check that this indeed is an equivalence relation onΩ D p (R d ).Perhaps more interestingly, we obtain the following result. For ease of notation we shall henceforth mean S(x) T when omitting the subscript from the signature. ) is injective up to tree-like equivalence, i.e., S(x) = S(y) iff x ∼ t y.Proposition A.1. Let p ≥ 1. The map S(•) : Ω D p (R d ) → G ( *Proof. The result follows from the continuous case upon realizing that S(x) = S(x) and analogously for y.In many cases, given a geometric càdlàg rough path x ∈ Ω D 0,p (R d ) with p ∈ [2, 3) and a path h ∈ D 1 ([0, T ], R e ) of bounded variation one is interested in constructing a new rough path y ∈ Ω D 0,p (R d+e ) such that the first level of y is given by y = (x, h). In the continuous case this can be done by using the level two information x 2 and dh ⊗ dh to fill in the corresponding terms in y 2 and using the well-defined Young cross-integrals to fill in the rest. The resulting level 2 rough path is called the Young pairing of x and h and we will denote it by y = P (x, h). The canonical example to keep in the mind is when h t = t, that is, we want to augment the rough path with an added time coordinate (see Def. A.8). In the càdlàg case one needs to be more careful in defining the appropriate Marcus lift. Definition A.7 (Definition 3.21 [7] ). Let x ∈ Ω D 0,p (R d ) with p ≥ 1 and h ∈ D 1 ([0, T ], R e ). Define the path z = (x, h) and the corresponding Marcus lift ẑ = (x, ĥ). The Young pairing of x and h is the p-rough path P (x, h) ∈ Ω D p (R d+e ) such that P (x, h) = P (x, ĥ) • η z where P (x, ĥ) is the usual Young pairing of a continuous rough path and a continuous bounded variation path (see Def. 9.27 in [17] ).We can then construct the time augmented rough path as the rough path obtained by the Young pairing with the simple continuous bounded variation path h t = t. It turns out that this pairing is continuous as a map fromΩ D 0,p (R d ) to Ω D 0,p (R d+1 ). Definition A.8. Let x ∈ Ω D 0,p (R d ).The time augmented version of x is the unique rough path x ∈ Ω D 0,p (R d+1 ) obtained by the Young pairing P (x, h) of x with the continuous bounded variation path h t = t. Proposition A.2. Let p ∈ [1, 3). Then, the map x → x is continuous and injective as a map fromΩ D 0,p (R d ) to Ω D 0,p (R d+1 ).Proof. Let X = Ω D 0,p (R d ) be a metric space when equipped with α p . Fix x ∈ X and let x n be a sequence of absolutely continuous paths converging in X to x. We shall first show that xn then converges to x. Since x n does not have any jumps and any reparameterisation of x n is still absolutely continuous, we may assume thatα p (x, x n ) = lim δ→0 d p (x δ , x n ) → 0for n → ∞. Define z = (x, h) and ẑd = (x δ , ĥδ ) the Marcus interpolation with η x,δ the reparameterisation such that z = ẑδ • η x,δ . Furthermore, let P (x, h) be the Young pairing of x and h. By definition, α p (P (x, h), P (x n , h)) = lim δ→0 σ p P (x δ , ĥδ ), P (x n , h)≤ lim δ→0 d p P (x δ , ĥδ ), P (x n , h) ≤ lim δ→0 C d p (x δ , x n ) + d 1 ( ĥδ , h) → 0for n → ∞ where C is just some generic constant depending only on p. The last inequality follows from 9.32 in [17] . Thus, if y ∈ X is such that α p (x, y) < ϵ, we can choose another sequence y n of absolutely continuous paths and N ≥ 1 large enough so that α p (P (x, h), P (y, h)) ≤ 2ϵ + α p (P (x n , h), P (y n , h)), α p (x n , y n ) ≤ 2ϵ for all n ≥ N . By Remark 3.6 in [7] , we then have that, up to choosing a large N , d p (x n , y n ) ≤ ϵ for all n ≥ N and therefore, once more appealing to Theorem 9.32 in [17] , α p (P (x n , h), P (y n , h)) ≤ 2Cϵ.In conclusion, α p (P (x, h), P (y, h)) ≤ 2(1 + C)ϵ which proves the result.Injectivity follows from [13] .The results of Section 3.3 hold in more generality. In fact, we can define Event RDEs similar to Definition 3.1 where the inter-event dynamics are given by Marcus RDEs driven by càdlàg rough paths. Utilizing the correspondence between solutions to Marcus RDEs and Marcus SDEs, it then follows that the results in the main body of the paper are a special case of the results given below. Definition A.9 (Event RDE). Let p ≥ 1 and N ∈ N be the number of events. Let x ∈ Ω D p (R d ) and f = (f 1 , . . . , f d ) be a family of Lip γ on R e with γ > p. Let E : R e → R and T : R e → R e be an event and transition function respectively. We say that y, (τ n ) N n=1 is a solution to the Event RDE parameterised by (y 0 , x, f, E, T , N ) if y T = y N T ,EQUATIONwith E(y n τn ) ̸ = 0 andEQUATIONEQUATIONExistence and uniqueness of solutions to Event RDEs is proven in the same way as for Event SDEs. Indeed, under the usual assumption that the vector fields f are Lip γ , for γ > p, a unique solution to (11) exists. In fact, the solution map y s × (s, t) → y t is a diffeomorphism for every fixed 0 ≤ s < t ≤ T (see, e.g., Theorem 3.13 in [7] ). It follows that we can iteratively define a unique sequence of solutions y n ∈ D p ([t n , T ], R d ). Finally, as mentioned in Remark A.2, if the driving rough path x is the Marcus lift of a semi-martingale, the inter-event solutions agree almost surely with the solutions to the corresponding Marcus SDE. Theorem A.2. Under Assumptions 3.1-3.2, there exists a unique solution (y, (τ n ) N n=1 ) to the Event RDE of Definition A.9. Furthermore, if x is the Marcus lift of a Brownian motion, the solution coincides almost surely with the solution to the corresponding Event SDE as given in Def. 3.1.Hence, the Event SDEs considered in the main text are special cases of Event RDEs driven by the Marcus lift of a Brownian motion. Yet, the more general formulation of Event RDEs allows to treat, using the same mathematical machinery of rough path theory a much larger family of driving noises such as fractional Brownian motion or even smooth controls. Also, since the driving rough path is allowed to be càdlàg, the model class given by Def. A.9 includes cases where the inter-event dynamics are given by Marcus SDEs driven by general semi-martingales.The proof of Theorem 3.2 presented below covers the case where (y, (τ n ) N n=1 ) is the solution to an Event RDE. Throughout we consider vector fields µ ∈ Lip 1 , σ ∈ Lip 2+ϵ and specialise to Event RDEs where the inter-event dynamics are given byEQUATIONwhere x ∈ Ω D p (R d ). The notation above deserves some clarification. One can define the vector field f = (µ, σ) and the Young pairing xt of x and h t = t. Assuming µ ∈ Lip 2+ϵ we can then view y n t as the unique solution to the Marcus RDEdy n t = f (y n t ) ⋄ xt .Alternatively, if one is not ready to impose the added regularity on the drift µ, one can view 13 as a RDE with drift as in Ch. 12 in [17] . To accommodate this more general case where the path driving the diffusion term might be 1) càdlàg and 2) is not restricted to be the rough path lift of a semi-martingale, we shall need the following two additional assumptions: Assumption B.1. For any n ∈ [N ], there exists a non-empty interval I n = (τ nδ n , τ n + δ n ) such that x is continuous over I n . In other words, the càdlàg rough path x, does not jump in small intervals around the event times (τ n ). Assumption B.2. For all 0 ≤ n ≤ N we define s n = τ nδ n /2 and t n = τ n+1 + δ n+1 /2. It holds that x ∈ Ω D 0,p ([s n , t n ], R d ), i.e., x is a geometric p-rough path on the intervals [s n , t n ]. Remark B.1. Note that Assumption B.1 trivially holds if x is continuous. Otherwise, it is enough to assume, e.g., that x is the Marcus lift of a finite activity Lèvy process. Furthermore, by the properties of the metric α p , if x is the canonical Marcus lift of a semi-martingale x ∈ D p ([s, t], R d-1 ), then there exists a sequence (x m ) of piece-wise linear pathsx m ∈ C 1 0 ([0, T ], R d-1 ) such that α p,[sn,tn] (x m , x) → 0 as m → ∞See, e.g. [7, Example 4.21] . The setting of Section 3.3 is therefore a special case of the setting considered here and Theorem 3.2 follows from the proof below. Proof. Recall that Φ(y, s, t; x) is the solution map or flow of the differential equationdy u = f (y u ) ⋄ dx u , y s = yevaluated at time t. The first step will be to prove continuity at y 0 . In particular, let y m 0 ∈ O approach y 0 for m going to infinity and denote the solutions to the corresponding Event RDEs by y m , (τ m n ) Nm n=1 . We claim that lim m→∞ N m = N andlim m→∞ τ m n = τ n , lim m→∞ y m,n-1 τ m n = y n τn .To see this, note that, by Theorem A.1, there exists a sequence λ m ∈ Λ of continuous reparameterisations such that |λ m | → 0 andsup (s,t)∈∆ T |Φ(y 0 , s, t; x) -Φ(y m 0 , λ m (s), λ m (t); x)| → 0 (14)for m → ∞. Note, furthermore, that Φ(y m 0 , s, t; x • λ m ) = Φ(y m 0 , λ m (s), λ m (t); x) for all (s, t) ∈ ∆ T . We let ỹm , (τ m n ) Nm n=1 be the solution to the Event RDE where (y 0 , x) is replaced by (y m 0 , x • λ m ). It suffices to prove that, for all 1 ≤ n ≤ N ,EQUATIONIndeed, since τ m n = λ -1 m (τ m n ) and |λ m | → 0, it then follows that τ m n → τ n for m → ∞. Furthermore, we have ỹm,n-1τ m n = y m,n-1 τ m n .We shall proof (15) using an inductive argument. We have thaty 0 t = Φ(y 0 , 0, t; x), ∀t ∈ [0, τ 1 ], ỹm,0 t = Φ(y m 0 , 0, t; x • λ m ), ∀t ∈ [0, τ m 1 ]. Now fix some 0 < ϵ < δ 1 where δ 1 is given in Assumption B.1. Note that |E(y 0 t )| > 0 for all t ∈ [0, τ 1ϵ] and therefore, by (14) , it follows that there exists an m 0 ∈ N such that, for all m ≥ m 0 ,inf t∈[0,τ1-ϵ] |E(Φ(y m 0 , 0, t; x • λ m ))| > 0 so that τ m 1 ≥ τ 1 -ϵ.Next, for some small 0 < η < ϵ, Assumption 3.4 and the Mean Value Theorem imply the existence of a+ η = r + η y 0 τ1 -(1 -r + η )y 0 τ1+η , and a - η = r - η y 0 τ1-η + (1 -r - η )y 0 τ1 with r + η , r - η ∈ (0, 1) such that E y 0 τ1+η = E y 0 τ1 + ∇E(a + η ) τ1+η τ1 µ(y 0 s )dy s , E y 0 τ1-η = E y 0 τ1 -∇E(a - η ) τ1 τ1-η µ(y 0 s )dy s ,But then, by Assumption 3.5 and the fact that E(y 0 τ1 ) = 0, for η small enough, E(y 0 τ1+η ) and E(y 0 τ1-η ) must lie on different sides of 0. Assumption B.1 and eq. ( 14) then yield the existence of am 1 ≥ m 0 such that τ m 1 ≤ τ 1 + η ≤ τ 1 + ϵ and inf t∈[0,τ1+η] |E(ỹ m,0 t )| > 0 for all m ≥ m 1 . It follows that τ m 1 → τ 1 . Finally, note that ỹm,0 τ m 1 -y 0 τ1 ≤ ỹm,0 τ m 1 -y 0 τ m 1 + y 0 τ m 1 -y 0 τ1 .Another application of (14) shows that the first term on the right hand side goes to 0 for m → ∞ and second term vanishes by Assumption B.1.To prove the inductive step, assume that (15) holds for i ≤ n. For all t ∈ [τ n , τ n+1 ] it holds that ỹm,n To prove part 3 we employ a similar induction argument to the one above. First, note that, by Theorem A.1, there exists a constant C > 0 not depending on x such that α p,[0,t0] y m,0 (a), y 0 (a) ≤ Cα p,[0,t0] (x m , x) .Since the latter term does not depend on y and goes to 0 for m going to infinity, we find that lim δ→0 τ m 1 (a) -τ δ 1 (a) + τ δ 1 (a) -τ 1 (a) . Now, let B 0 = B r (y 0 ). Since τ 1 (a) ∈ (τ 1 -δ 1 /2, τ 1 + δ 1 /2) for all a ∈ B 0 and x is continuous over this interval, it follows that τ δ 1 (a)τ 1 (a) goes to 0 as δ → 0 for each a ∈ B 0 . Furthermore, by definition of the metric α p , eq. ( 16), and the fact that y m,0 0 (a) = a = y 0 0 (a), for each a ∈ B 0 , a similar argument as the one employed in the beginning of the proof then shows that |τ τ δ 1 (a) (a)y 0 τ1(a) (a) and taking the limit as δ → 0 and then the supremum over x ∈ B 0 on both sides, we can argue in exactly the same way to show that part 3 holds for n = 1. We can then argue by induction, just as in the first part of the proof, to show that it holds for all subsequent event times as well. Thus, the set B 0 satisfies all the stated requirements.Lemma B.2. Let Assumption B.1 hold and x m be as in Assumption B.2. Then, for all n ∈ [N ] and p ′ > p, limm→∞ d p ′ ,[s,t] (x m , x) = 0, for any τ n -δ n /2 ≤ s < t ≤ τ n + δ n /2.Proof. Fix some n ∈ [N ], p ′ > p and τ nδ n /2 ≤ s < t ≤ τ n + δ n /2. Note that, for any continuous reparameterisation λ ∈ Λ, m ∈ N, and δ > 0, it holds thatd p ′ ,[s,t] (x m , x) ≤ d p ′ ,[s,t] (x m , x m • λ) + d p ′ ,[sn,tn] (x m • λ, xδ ) + d p ′ ,[s,t] (x δ , x),where xδ is the Marcus interpolation of x over the interval [s n , t n ]. Taking the infimum over λ ∈ Λ and the limit as δ → 0 on both sides, we obtaind p ′ ,[s,t] (x m , x) ≤ α p ′ ,[sn,tn] (x m , x) + lim δ→0 d p ′ ,[s,t] (x δ , x).The first term on the right hand side goes to 0 as m → ∞ by Assumption B.2. Furthermore, since, byAssumption B.1, x is continuous on (τ n -δ n , τ n + δ n ), it follows that d ∞,[s,t] (x δ , x) goes to 0 for δ → ∞.But the result then follows from Proposition 8.15 and Lemma 8.16 in [17] .Proof of Theorem 3.2.Step 1:Assume that x ∈ C 1 ([0, T ], R d-1 ). By [17, Theorem 4.4] , the Jacobian ∂y 0 t exists and satisfies (7) for all t ∈ [0, τ 1 ). We shall prove that relations ( 6) and ( 7) hold for all n ∈ [N ] by induction. Thus, assume that ∂y k t and ∂τ k exist for all t ∈ [τ k , τ k+1 ) and k ≤ n -1 and satisfy the stated relations. To emphasise the dependence on the initial condition, we will sometimes use the notation y n = y n (y 0 ) and τ n = τ n (y 0 ) for the solution of the Event RDE started at y 0 . We want to show that, for arbitrary h ∈ R e , the following limitslim ϵ→0 τ ϵ n -τ n ϵ and lim ϵ→0 y n,ϵ t -y n t ϵ for t ∈ [τ n , τ n+1 )exist and satisfy the stated expressions, where τ ϵ n = τ n (y 0 + hϵ) and y n,ϵ = y n (y 0 + hϵ). For any ϵ > 0, because E is continuously differentiable, the Mean Value Theorem implies that there exists c ϵ ∈ R e on the line connecting y n-1 τn to y n-1 τ ϵ n and another c ′ ϵ ∈ R e on the line connectingy n-1,ϵ τ ϵ n to y n-1 τ ϵ n such that E y n-1 τn = E y n-1 τ ϵ n + ∇E(c ϵ ) y n-1 τn -y n-1 τ ϵ n = E y n-1 τ ϵ n + ∇E (c ϵ ) µ(y n-1 τn )(τ n -τ ϵ n ) + σ y n-1 τn (x τn -x τ ϵ n ) + o(|τ n -τ ϵ n |) , E y n-1,ϵ τ ϵ n = E y n-1 τ ϵ n + ∇E(c ′ ϵ ) y n-1,ϵ τ ϵ n -y n-1 τ ϵ n = E y n-1 τ ϵ n + ∇E(c ′ ϵ ) ϵ ∂y n-1 τn h + o(ϵ) ,where the last equality follows from the induction hypothesis. We haveE(y n-1 τn ) = 0 = E(y n-1,ϵ τ ϵ n). Thus, by rearranging, we find thatτ ϵ n -τ n ϵ = - ∇E(y n-1 τn )∂y n-1 τn h ∇E(y n-1 τn ) µ(y n-1 τn ) + σ(y n-1 τn ) xτ n -x τ ϵ n τn-τ ϵ n + o(1) = - ∇E(y n-1 τn )∂y n-1 τn h ∇E(y n-1 τn )µ(y n-1 τn ) + o(1)where the second equality follows from Assumptions 3.4 and 3.5.Assume for now that τ ϵ n < τ n . By another application of the Mean Value Theorem, there exists c ϵ ∈ R e on the line connecting y n-1 τn to y n-1 τ ϵ n such thaty n,ϵ τn -y n τn = y n,ϵ τn -T y n-1 τn = y n,ϵ τn -T y n-1 τ ϵ n -∇T (c ϵ )(y n-1 τn -y n-1 τ ϵ n ) = y n,ϵ τ ϵ n + µ(y n,ϵ τn )(τ n -τ ϵ n ) + σ y n,ϵ τn (x τn -x τ ϵ n ) -T y n-1 τ ϵ n -∇T (c ϵ ) µ(y n-1 τn )(τ n -τ ϵ n ) + σ y n-1 τn (x τn -x τ ϵ n ) + o(|τ n -τ ϵ n |) = T y n-1,ϵ τ ϵ n -T y n-1 τ ϵ n + µ(y n,ϵ τn ) -∇T (c ϵ )µ(y n-1 τn ) (τ n -τ ϵ n ) + σ(y n τn ) -∇T (c ϵ )σ(y n-1 τn ) (x τn -x τ ϵ n ) + o(|τ n -τ ϵ n |) Therefore y n,ϵ τn -y n τn ϵ = ∇T (y n-1 τn )∂y n-1 τn h + µ(y n τn ) -∇T (y n-1 τn )µ(y n-1 τn ) ∂τ n h + o(1)where we used Assumption 3.3, the chain rule and the existence of ∂τ n . Finally, for any t ∈ (τ n , τ n+1 ], equation ( 7) follows from the fact that we can write y n t = Φ(y n s , s, t, x) for all τ n ≤ s < t. In particular, by the chain rule, we find that∂y n t = ∂ y n s Φ(y n s , s, t)∂y n s s=τn = ∂ y n τn y n t [∂y n s ] s=τn .Step 2: Consider now the general case of x ∈ Ω p (R e ) and let (y m , (τ m nm ) Nm nm ) denote the solution to the Event RDE where x is replaced by the piece-wise linear approximation x m . With ∂y n,m t and ∂τ m n denoting the corresponding derivatives, we saw in the previous step that both exist and satisfy ( 6)- (7) . We let R n t and ρ n denote the right hand side of ( 7) and ( 6) respectively. This step consists of proving that, for n ∈ [N ] and t ∈ (τ n , t n ),EQUATIONand for some open ball B 0 around y 0EQUATIONBy Lemma B.1 and continuity of T we have that T (y m,n-1 τ m n ) converges to T (y n-1 τn ) and τ m n converges to τ n as m → +∞. Then, because y n t = Φ(T (y n-1 τn ), τ n , t; x) and y m,n t = Φ(T (y m,n-1EQUATIONfollows from, Lemma B.2 and Corollary 11.16 in [17] . In fact, since B 0 was constructed in Lemma B.1 in such a way that τ n (a) < t n for all a ∈ B 0 we also get thatlim m→∞ sup a∈B0 ∂ y n τn (a) (a) Φ y n τn(a) (a), τ n (a), t; x -∂ y m,n τ m n (a) (a) Φ y m,n τ m n (a) (a), τ m n (a), t; x m = 0.by the same corollary in [17] . Thus, to prove (18) , it suffices to show that, for all n ∈ {1, ..., N },lim m→∞ sup a∈B0 ∂y m,n-1 τ m n (a) (a) -R n-1 τn(a) (a) = 0.We shall prove it using another inductive argument starting with n = 1. In this case it suffices to show that limm→∞ sup a∈B0 ∥∂ a Φ (a, 0, τ 1 (a); x) -∂ a Φ (a, 0, τ m 1 (a); x m )∥ = 0.By [7, Theorem 3.3] we know that the above holds if τ m 1 (a) and τ 1 (a) are replaced by τ 1 + δ 1 /2. Now let Φ -1 be the reverse of the flow map Φ, that is,Φ -1 (a 1 , s, t; x) = a 0 ⇔ Φ(a 0 , s, t; x) = a 1 .From Lemma B.1 it follows that y 0 τ1(a) (a) = Φ -1 (y 0 t0 (a), τ 1 (a), t 0 ; x) and, for m large enough, y m,0 τ m 1 (a) (a) = Φ -1 (y m,0 t0 (a), τ m 1 (a), t 0 ; x m ). But the result then follows from Lemma B.2 and [17, Corollary 11.16 ]. To prove the inductive step, assume that (18) holds for all i ≤ n -1. Again, by inspecting ( 6) and ( 7) and using the inductive assumption, one finds that it is enough to show thatlim m→∞ sup a∈B0 ∂ y n-1 τ n-1 Φ y n-1 τn-1 , τ n-1 , τ n ; x -∂ y m,n-1 τ m n-1 Φ y m,n-1 τ m n-1 , τ m n-1 , τ m n ; x m = 0,where we suppressed the dependence on a for notational simplicity. This is done exactly as for y 0 and completes the proof of Step 2.Step 3: The third and final step is to combine Step 1 and 2 to finish the proof. So far we have proven that 1) the theorem holds for continuous paths of bounded variation and 2) (τ m n , y m,n t ) converges to (τ n , y n t ) and (∂τ m n (a), ∂y m,n t (a)) converges uniformly to (ρ n (a), R n t (a)) over a ∈ B 0 for all t ∈ (τ n , t n ) and n ∈ [N ]. From these results it immediately follows that (τ n (a), y n t (a)) is differentiable at a = y 0 with derivatives given by (ρ n (y 0 ), R t (y 0 )) for all t ∈ (τ n , t n ). What is left to show then, is that this also holds for all other t. But this follows immediately from the chain rule upon realizing that, for any τ n < s < τ n + δ n /2 < t < τ n+1 ,y n t = Φ(y n s , s, t; x) ⇒ ∂y n t = ∂ y n s Φ(y n s , s, t; x)R n s (y 0 ) = R n t (y 0 ).We give here a brief outline of some of the most central concepts related to kernel methods. For a more in-depth introduction we refer the reader to [45, 55, 3] . Let X be a topological space. We shall in this paper only be concerned with positive definite kernels, that is, symmetric functions k : X × X → R for which the Gram matrix is positive definite. To such a kernel one may associate a feature map X → R X such that x → k x = k(x, •). A reproducing kernel Hilbert space (RKHS) is a Hilbert space H ⊂ R X such that the evaluation functionals, ev x : f → f (x), are bounded for each x ∈ X . For all positive definite kernels there is a unique RKHS H ⊂ R X such that f (x) = ⟨k x , f ⟩ H for all f ∈ H and x ∈ X . This is also known as the reproducing property. Furthermore, with H denoting the linear span of {k x | x ∈ X }, it holds that H = H, i.e., H is dense in H. Two important properties of kernels are characteristicness and universality. Definition C.1. Let k : X × X → R be a positive definite kernel. Denote by H the linear span of {k x | x ∈ X } and let F ⊂ R X be a topological vector space containing H and such that the inclusion map ι : H → F is continuous.• We say that k is universal to F if the embedding of ι : H → F is dense.• We say that k is characteristic to F ′ if the embedding µ :F ′ → H ′ , D → D| H is injective Remark C.1.This definition is the one used in [8] and is more general then the one usually encountered. Note that in many cases (all the cases considered here, in fact) F ′ will contain the set of probability measures on X in which case k being characteristic implies that the kernel mean embedding µ→ E X∼µ k X (•) is injective. Remark C.2.Often times, instead of starting with the kernel function k and then obtaining the RKHS, one starts with a feature map F : X → H into a RKHS and then defines the kernel as the inner product in that Hilbert space, i.e., k(x, y) = ⟨F (x), F (y)⟩ H . In such cases, it makes sense to ask whether there are equivalent notions of F being universal and characteristic. This is indeed the case and the definition is almost the same as above. We refer to Definition 6 in [8] for a precise statement.The definition of the signature kernel requires an initial algebraic setup. Let ⟨•, •⟩ 1 be the Euclidean inner product on R d . Denote by ⊗ the standard outer product of vector spaces. For any n ∈ N, we denote by ⟨•, •⟩ n on (R d ) ⊗n the canonical Hilbert-Schmidt inner product defined for any a = (a 1 , . . . , a n ) and b= (b 1 , . . . , b n ) in (R d ) ⊗n as ⟨a, b⟩ n = n i=1 ⟨a i , b i ⟩ 1 . The inner product ⟨•, •⟩ n on (R d ) ⊗n can then be extended by linearity to an inner product ⟨•, •⟩ on T ((R d )) defined for any a = (1, a 1 , . . .) and b = (1, b 1 , . . .) in T ((R d )) as ⟨a, b⟩ = 1 + ∞ n=1 ⟨a n , b n ⟩ n . To begin with, let X = D 1 ([0, T ], R d ).If x ∈ X is càdlàg path, we can define the Marcus signature in the spirit of Marcus SDEs [42, 43] as the signature of the Marcus interpolation of x. This interpolation, denoted by x, is the continuous path on [0, T ] obtained from x by linearly traversing the jumps of x over added fictitious time r > 0 and then reparameterising so that the path runs over [0, T ] instead of [0, T + r]. The general construction is given in Appendix A. If x is continuous, x and x coincide; thus, without any ambiguity, we can define the Marcus signature S(x) of a general bounded variation càdlàg path as the tensor series described above, but replacing x with x (see also the definition in A.2).Since the signature is invariant to certain reparameterisations (Proposition A.1), it is not an injective map. Injectivity is a crucial property required to ensure characteristicness of the resulting signature kernel that we will introduce next. One way of overcome this issue is to to augment a path x with a time coordinate resulting in the path x = (x, t) 5 . The Marcus signature kernel is then naturally defined as the map k : X × X → R such that k(x, y) = ⟨S(x), S(ỹ)⟩ for any x, y ∈ X . As stated in Theorem C.1, this kernel is universal on compact subsets K ⊂ X and, equivalently, characteristic to the space of regular Borel measures on K. However, these properties do not generalize to the whole space C b (X , R) of bounded continuous functions from X to R.In [8] the authors address this issue in the case of continuous paths by introducing the so-called robust signature. They define a tensor normalization as a continuous injective mapΛ : T R d → a ∈ T R d | ∥a∥ ≤ Rfor some R > 0 and such that Λ(a) = (a 0 , λ(a)a 1 , λ(a) 2 a 2 , . . . ) for some λ : T R d → (0, ∞). Now, let p ∈ [1, 3) and take C 1 0 (R d ) to be the space of absolutely continuous functions onR d . Recall that Ω D 0,p (R d ) is the closure of C 1 0 (R d ) in Ω D p (R d) under the metric α p . Throughout we let X = Ω D 0,p (R d ) be a metric space equipped with α p . Naturally, we can then define the signature kernel on X by k(x, y) = ⟨S(x), S(ỹ)⟩ and, similarly, the robust signature kernel k Λ (x, y) = ⟨Λ(S(x)), Λ(S(ỹ))⟩ where Λ is a tensor normalisation.Theorem C.1. Let p ≥ 1, Λ a tensor normalization, and K ⊂ X compact under α p . Then, (i) The signature kernel k is universal to F = C(K, R) equipped with the uniform topology and characteristic to the dual F ′ , the space of regular Borel measures on K.(ii) The robust signature kernel k Λ is universal to F = C b (X , R) equipped with the strict topology and characteristic to the dual F ′ , the space of all finite Borel measures on X .Proof of Theorem C.1. Part (i) follows directly from the proof of Proposition 3.6 in [13] . For part (ii) we shall proof that the feature map F = Λ • S is universal and characteristic. The result then follows from Proposition 29 in [8] . We start by defining P = X / ∼ t where the equivalence relation ∼ t is defined in Appendix A.2. We equip P with the topology induced by the embedding S : P → T ((R d+1 )). By Proposition A.1, F is a continuous and injective map from P into a bounded subset of T ((R d+1 )). Thus,H = {⟨ℓ, F ⟩ | ℓ ∈ T ((R d+1)) ′ } is a subset of F that separates points. Furthermore, since F takes values in the set of group-like elements, H is a subalgebra of F (under the shuffle product). It then follows from Theorem 7 and Theorem 9 in [8] that F is universal and characteristic. The fact that F ′ is the space of all finite Borel measures on X is part (iii) of Theorem 9 in the same paper. Finally, as per Appendix A.3, the map x → x is a continuous and injective embedding of X into P from which the result then follows.With d k denoting the MMD for a given kernel k : X × X → R, the following is a direct consequence of Theorem C.1.Corollary C.1. Let p ≥ 1, Λ a tensor normalization, and K ⊂ X compact under α p . Then, d k is a metric on M(K) and d kΛ is a metric on M(X ).In the general SSNN model, Theorem 3.2 gives the following result. Proposition D.1. Fix some weight w ij ∈ w, a neuron k ∈ [K] and let G k t denote the gradient of (v k , i k ) wrt. w ij at time t. Furthermore, define γ : {0, 1} → R 2 such that γ 0 = (µ 1 , -µ 2 )w lk , γ 1 = (µ 1 , 0)v reset , and let Γ ∈ R 2×2 be the drift matrix in the inter-spike SDE of (v k , i k ). Then,EQUATION)where e n ∈ R 2 is the n'th unit vector, l is the neuron in Pa k ∪ {k} with the most recent spike time before t, and we denote this spike time by s. If t is a spike time of neuron k it therefore follows thatEQUATION)where t prev is the previous spike time of neuron k. In the case of a deterministic SNN, formula (20) is replaced byEQUATION)Proof. Throughout we fix some t > 0 and let s < t denote most recent event time preceding t with l the index of the neuron firing at time s. We define the process dw t = 0dt with w 0 = w ij and with a slight abuse of notation we shall write y k t = (v k t , i k t , s k t , w t ). We will leave out the event index n for notational simplicity. Since y k t depends on y s only through y k s and ∇T l is block diagonal, a direct consequence of eq. ( 7) isG k t = (I 0) ∂ y k s y k t ∇T l j (y k s-)∂ wij y k s--µ(y k s ) -∇T k l (y k s-)µ(y k s-) ∂ wij s where µ(v, i, s, w) = (µ 1 (i -v), -µ 2 i, λ(v), 0). If l ∈ Pa k ∪ {k}, then T k l = id and therefore G k t = (I 0) ∂ y k s y k t ∂ wij y k s .One can then reapply the formula above until l ∈ Pa k ∪ {k}. By the flow property, it follows that we may assume without loss of generality that l ∈ Pa k ∪ {k}. This leaves us with two cases. We definez k t = (v k t , i k t ) so that (I 0) ∂ y k s y k t = ∂ z k s z k t and ∂ wij z k t = G k t . Furthermore, let a = δ il δ jk .Case 1, l ∈ Pa k : In this case T k l (v, i, s, w) = (v, i + aw + (1a)c, s, w) where c is a constant. As a result∂ z k s z k t ∇T k l (y k s-)∂y k s-= ∂ z k s z k t G k t + a∂ i k s z k t , ∂ z k s z k t µ(y k s ) -∇T k l (y k s-)µ(y k s-) = ∂ z k s z k t γ 0 , In total, G k t = ∂ z k s z k t G k t -γ 0 ∂ wij s + ae 2 . Case 2, l = k: In this case T k l (v, i, s, w) = (v -v reset , i, log u -α, w) so that ∂ z k s z k t ∇T k l (y k s-)∂y k s-= ∂ z k s z k t G k t , ∂ z k s z k t µ(y k s ) -∇T k l (y k s-)µ(y k s-) = ∂ z k s z k t γ 1 , and, thus, G k t = ∂ z k s z k t G k t -∂ z k s z k t γ 0 ∂ wij s.Note that z k t is an Ornstein-Uhlenbeck process initialized at z k s and with drift and diffusion matricesΓ = -µ 1 µ 1 0 -µ 2 , Σ = σ 1 0 0 σ 2 .As a result, we can directly compute ∂ z k s z k t = e (t-s)Γ . This proves that eq. ( 19) holds. Eq. ( 20) then follows directly from ( 6) and the fact that E k (y) = s k .From this the results of Section 4.4 follow since the terms ∂ wij s vanish whenever s is the spike time of a neuron l that is not a descendant of neuron j. Thus, equation (19) only includes terms depending on the activity of the pre-and post-synaptic neuron. In particular, there is no need to store the gradient path G k t for each combination of neuron k and synapse ij, but each neuron only needs to keep track of the paths for its incoming synapses. This reduces the memory requirements from the order of K 3 to only K 2 (which is needed anyway to store the weight matrix). In general, the gradient paths can be approximated by simply omitting the terms ∂ wij s.E.1 Input current estimation 6 For each combination of sample size and σ we sample a data set of spike trains using Algorithm 1 with N = 3, i.e., up until the first three spikes are generated. We use diffrax to solve the inter-Event SDE with a step size of 0.01 and the numerical solver is the simple Euler-Maruyama method. We then sample an initial guess c ∼ Unif([0.5, 2.5]) and run stochastic gradient descent using the approach described in 4.1. That is, for each step, we generate a batch of the same size as the sample size and use d k to compare the generated batch to the data. For each step we also compare the absolute error between the average spike time of the first three spikes of the generated sample to a hold a out test set of the same size as the sample. We use the RMSProp algorithm with a decay rate of 0.7 and a momentum of 0.3 which we found to work well in practice. The learning rate is 0.001. The experiment was run locally on CPU with an Apple M1 Pro chip with 8 cores and 32 GB of ram. The entire experiment took approximately 3-6 hours to run. For the exact details of this experiment we refer to the notebook snnax/notebooks/single_neuron.ipynb in the supplementary material.As above, for each sample size D ∈ {256, 512, 1024} we sample a data set of spike trains using Algorithm 1 with T = 1 and with the same differential equation solver setup as above. Thus, in this case, the number of spikes varies across each sample path. The parameters are chosen as follows:• v reset = 1.2 • λ(v) = exp(5(v -1) • µ = (6, 5) • σ = I 2 /4For each sample size the data was generated using the same randomly sampled weight matrix w which represents a feed-forward network of the dimensions described in Section 4 and which was constructed as follows: for the weight matrix from layer l to layer l + 1, say w l , we sample each entry from Unif([0.5, 1.5]) and then normalize by 3/K l where K l is the number of neurons in layer l. The normalisation makes sure that the spike rate for the neurons in each layer is appropriate.For each data set (each sample size) we then train a spiking neural net of the same network structure to match the observed spike trains. This is done using stochastic gradient descent with a batch size of B = 128 and by computing d k on a generated batch and a batch sampled from the data set at each step. In order to avoid local minimums 7 we match the number of spikes between the generated spike trains and the ones sampled from the data set. Also, we sample from the data set without replacement so that we loop through the whole data set every D/B steps. We run RMSProp for 1500 steps with a momentum of 0.3 and a learning rate of 0.003 for the first 1000 steps and 0.001 for the last 500 steps. This experiment was run in the cloud using Azure AI Machine Learning Studio on a NVIDIA Tesla V100 GPU with 6 cores and 112 GB of RAM. The entire experiment took around 12-16 hours to run. For the exact details we refer to the notebook snnax/notebooks/spiking_neural_net.ipynb in the supplementary material. Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Justification: There is a dedicated section (Section 5) clearly discussing what is lacking and what can be explored in future work.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.The point here is actually a little more subtle. It is in fact possible to obtain exact gradients using the adjoint method as long as one uses reversible (or adjoint) solvers in the forward pass.These are standard regularity conditions to ensure existence and uniqueness of a strong solution.For simplicity we only compute an approximation of the true MMD by truncating the signatures at depth 3 and taking the average across the batch/sample size.Note that in the case of stochastic SNNs this term is not necessarily well-defined since semi-martingales are in general not differentiable wrt. time.If x is a càdlàg rough path, this is done via a Young pairing which results in a càdlàg p-rough path, x, where the first level is given by (xt, t). For more information, we refer to Appendix A.3.For an updated version of some of the experiments and a general implementation of SSNNs in diffrax we refer to github.com/cholberg/snnax.Note that the loss landscape is inherently discontinuous since whenever the parameters are altered in such a way that an additional spike appears, the expected signature will jump. 


Decision: This paper describes a novel formalism for computing gradients within spiking neural networks using stochastic differential equations.  All four reviewers felt that it was above the threshold for acceptance, and I'm pleased to report that it has been accepted to NeurIPS.  Congratulations!  Please revise the manuscript according to the reviewer comments and discussion points

Reviews: review: paper_summary: paper_summary: This paper introduces a general mathematical framework to model stochastic spiking neural networks (SSNN) as stochastic differential equations with event discontinuities, and identifies sufficient conditions ensuring the existence of gradients. With a newly defined loss function, SSNNs can be trained as generative models with an end-to-end autodifferentiable solver. Some empirical verifications demonstrate the effectiveness of the proposed method, and also, discussions on the connection with bio-plausible learning algorithms are provided.
strength_weakness: strength_weakness: 1. There are several assumptions in the analysis. It would be better to have more discussions on whether these assumptions (e.g., Assumptions 3.1 and 3.2) can hold in common conditions.

2. This paper only conducts experiments on two toy settings. It would be better to discuss how to apply the method to more application conditions.
questions: questions: What is the computational complexity of SDESolveStep and RootFind?
limitations: limitations: The authors discussed limitations in Section 5.
paper_summary: paper_summary: The paper introduces a mathematical framework using rough path theory to model stochastic spiking neural networks as stochastic differential equations with event discontinuities, driven by càdlàg rough paths.  This framework accommodates potential jumps in both solution trajectories and driving noise. Furthermore, the authors introduce Marcus signature kernels to extend continuous signature kernels to càdlàg rough paths and using them to define a general-purpose loss function for training stochastic SNNs as generative models.
strength_weakness: strength_weakness: I am not an expert in rough path theory and I don't see any major weaknesses.
questions: questions: 1. How does the performance of the proposed gradient-based training compare to other established methods in the literature, such as those in [1,2,3]? More exhaustive empirical evaluations would be valuable and this would help to understand the benefits and limitations of this approach in practice.

2. Can the framework be extended to the case of single spike scenario (for instance in [1,2,3])? 

[1] Göltz, J., Kriener, L., Baumbach, A. et al. Fast and energy-efficient neuromorphic deep learning with first-spike times. Nat Mach Intell 3, 823–835 (2021).

[2] I. M. Comsa, K. Potempa, L. Versari, T. Fischbacher, A. Gesmundo and J. Alakuijala, "Temporal Coding in Spiking Neural Networks with Alpha Synaptic Function," ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Barcelona, Spain, 2020, pp. 8529-8533, doi: 10.1109/ICASSP40776.2020.9053856.

[3] Stanojević, Ana et al. “An Exact Mapping From ReLU Networks to Spiking Neural Networks.” Neural networks : the official journal of the International Neural Network Society 168 (2022): 74-88.
limitations: limitations: There is no potential negative societal impact.
paper_summary: paper_summary: This work develops a mathematical framwork to compute the gradients of stochastic spiking neural networks, or more generally Event SDEs. The proposed framework is an alternative to existing surrogate gradient frameworks and an extension of prior adjoint based work (e.g. EventProp) in the presence of stochasticity. Jax-based code is provided with the submission.
strength_weakness: strength_weakness: - Use of jargon and mathematical assumptions which are sometimes not sufficiently connected with the studied problem. This is a missed opportunity IMO to connect with the relevant communities. e.g. I have little clue what this means:
  "train NSDEs non-adversarially using a class of maximum mean discrepancies (MMD) endowed with signature kernels [51] indexed on spaces of continuous paths as discriminators". 
 Why adversarial? what is a continuous signature kernel? 
 Maybe I'm the wrong reader for this paper, but then this would likely apply to most readers interested in training SNNs.

- Lack of a benchmark and experimental comparison to prior work based on eventprop or surrogate gradients. The framework is written in jax, and should in principle allow a scalable simulation, why not run on a simple SNN benchmark? The provided experiments are overly simple and leaving this open as a vague limitation statement is not enough.
questions: questions: - Please explain some ambiguous sentences and less known facts:
  - L211: "If the diffusion term is absent we recover the gradients in [6]" from which expression to we recover the gradients in [6]?
  - Give at least a one sentence description of cadlag (right continuous, left limit). Why is it relevant to the modeling of SNNs?
  - What is a continuous signature kernel? Can you give an intuition without having to read another technical paper?
- Definition 3.1and its eq (3-5) are to my understanding a general version of eq (2). Is it really necessary to define these two models separately? I would recommend using the same notation if possible, or to the very least to connect the expressions across the two definitions.
- Can you solve at least one practical SNN training problem?
- How does the training effort (wall time, memory) compare to existing methods? Is this an impediment to test on a benchmark?
limitations: limitations: Some limitations are discussed in the conclusions, but leaving a simple SNN benchmark to future work is not a limitation, but a shortcoming of the manuscript in my opinion.
paper_summary: paper_summary: Spiking neural networks face the problem of non-differentiability of loss function due to the Heaviside activation function, which works as the spiking function. To back-propagate the loss through the network, Heaviside is replaced with surrogate functions as a workaround. 

The paper establishes the theoretical existence of the exact gradient for stochastic spiking neural networks, where stochasticity is introduced in spike timing, and the gradient of spike timing is computed with respect to the model parameters. It uses the existing theories of stochastic differential equations (rough path theory) to handle discontinuous spiking events and introduces a novel kernel (Marcus signature kernel) to handle such discontinuities. They further provide an algorithm that can perform parameter estimation of a two-layer network on toy examples.
strength_weakness: strength_weakness: The paper demonstrates it implementation through a toy parameter estimation problem on an one-hidden layer feed forward network of dimension $4 \times 16 \times 2$, where the network parameters are estimated. Such a network, having a tiny input dimension, is not readily applicable to realistic datasets. Further, rather than implementing the algorithm as a parameter search (as done in standard SNN training), it is implemented as a parameter estimation problem. The algorithmic aspects of the paper are somewhat open; the paper discusses several possible avenues, such as "discretize-then-optimize" vs. "optimize-then-discretize" frameworks. The authors must discuss the computation complexity of the proposed EventSDESolve algorithm, which uses the discretize-then-optimize approach.
questions: questions: **Doubt:** Assumptions 3.3, 3.4 and 3.5 require the event and transition function to be continuously differentiable [line 194]. How do such functions model discontinuous spikes? In contrast, line 229 mentions that trajectories and spike times are almost surely differentiable. Please clarify.

**Measure of Randomness:** It is shown that the zeroth-order derivative of the Heaviside function approximates a surrogate function depending upon the noise distribution [1]. How, in the present case, does the variance of the introduced noise (in the Brownian motion) affect the exact gradient? Does it recover the gradient of Heaviside in a limiting case? Which of the two noises (Brownian motion vs. transition noise [line 166]) used in the SLIF modelling is essential to obtain the exact gradient -- are either of them a mathematical convenience? 

**Computational Feasibility:** Please discuss the computational complexity of EventSDESolve. How does the computational complexity of the gradient scale with the number of parameters? Can the algorithm handle real datasets? Why did the authors choose to solve a parameter estimation problem (through estimation of input current) instead of standard training of the network parameters? 

[1] Mukhoty, Bhaskar, et al. "Direct training of snn using local zeroth order method." Advances in Neural Information Processing Systems 36 (2024).
limitations: limitations: Yes, authors do discuss some of the limitations of their work.

Comments: comment: Thank you for your comments and rebuttal. It would be valuable to include the points in the general response from the authors in the paper for additional context. I will keep my score.
I would like to thank the authors for their detailed responses and discussions, and I keep my score.
I thank the authors for the clarification presented. It could be beneficial for future readers if some of the concerns raised by the reviewers and the corresponding responses could be summarized in the paper's appendix.
Thanks you for the clarifications, I maintain my score.