Exact Gradients for Stochastic Spiking Neural Networks Driven by Rough Signals
Authors: Christian Holberg
Cristopher Salvi
Venue: NeurIPS 2024 poster (2024)
Keywords: stochastic differential equations
spiking neural networks
backpropagation
rough paths
signatures
Primary area: neuroscience_and_cognitive_science

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.

Exact Gradients for Stochastic Spiking Neural Networks Driven by Rough Signals
Christian Holberg Department of Mathematics University of Copenhagen hoholba holberg@math ku dk
Cristopher Salvi Department of Mathematics Imperial College London sasalvi @imperial allall.. ac uk
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 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 introduce a new class of signature kernels
how these gradients satisfy recursive relation. Furthermore, general-purpose loss function defined by means 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 part of the diffrax library. Our framework is, our knowledge, the first enabling gradient- based training of SSNNs with noise affecting both the spike timing and the network's 's dynamics.
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-frr (SLIF) model being among the most popular choices (6] In its simplest form, given some continuous input current ig on [0, T). the dynamics of a single SLIF process describing the membrane potential well
neuron consist of an Crnstein-Uhlenbed threshold for spike triggering and a resetting mechanism [33] In particular, between spikes, the dynamics of the membrane potential is given by the following SDE
du; = p (i+ v2) dt + adB,,
(1)
where μ > 0 is parameter and B, is a standard Brownian motion. The neuron spikes whenever the membrane potential hits the threshold & upon which is reset 0. Alternatively, one can model the spike times as Poisson process with intensity X IR + R. depending on the membrane potential Ug. A common choice X(v) cxp((r 4)/B) 50
24
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 's parameters using classical calculus rules. This issue exacerbated when the dynamics implicit function theorem, used for instance in 6 is no longer valid.
stochastic in which case the usual argument relying the
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

1.1 Contributions
In this paper, we introduce 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 that of rough path theory (38] a modern branch
of stochastic analysis providing robust solution theory for stochastic dynamical systems driven by noisy, possibly discontinuous, rough signals. Although Brownian motion a prototypical example,
these signals can be far more irregular (or rougher) than semimartingales [7 16
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 '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 general-purpon loss function enabling the training of SSNNs as generative models. We provide an end-to-end autodifferentiable solver for Event SDEs (Algorithm D 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 numerical solver (not necessarily based on a Euler scheme). In fact, previous solutions are based either on surrogate gradients [46] or follow an optimise- then- don--iicretise approach deriving adjoint equations [56), the latter yielding exact gradients only the scenario where solutions are available in closed form and not approximated via numerical solver T Finally, we discuss how our results lead to bioplausible learning algorithms akin to e-prop
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: fo )t started at yo IR using parametric Lipschitz vector field fo R - R2, usually given by a neural network 5) Similarly, a neural stochastic differential equation (NSDE) an SDE of the form dy, = poy )dl + +o(yy))d driven by a d-dimensional Brownian motion B. started at yn E IR and with parametric vector field p4 : IR + Re and og R Rexd that are Lip² and Lip2+e continuous respectively 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) [4 22). For an account on applications of rough path theory to machine learning see 14 [1].
Training techniques for NSDES Training 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 W-asserstein distance to train NSDE as a Wasserstein-GAN []]
Latent SDEs (36] train NSDE with respect to the KI. divergence via variational inference and can be interpreted as variational autoencoders.
23 the authors propose to train NSDEs aov-avversarially using 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 48 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.
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
Backpropagation through NSDEs
'The point here is actually a little more subtle. It in fact possible to obtain exact gradients using the adjoint method as long as one uses reversible (or adjoint) solvers the forward pass.
"hhee are standard regularity conditions to ensure existence and uniqueness of strong solution.

---

first option is simply to backpropagate through the solver's 's internal operations. This method is known as wwrreetsee-tee-ppttniii it is generally speaking fast evaluate and produces accurate gradients. every internal operation of the solver must be recorded. A second computes gradients by deriving backwards-in-time
but it is memory-ineffiiiett, approach, known as ppiiiis--teen-dieree 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 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 When the differential equation is an ODE, there is rich literature on sensitivity analysis aimed at computing derivatives using the implicit function theorem
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 [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 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. where the backward pass uses custom gradients for the non 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 -heee-scete-iiseeeeeeeeeee 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 dittediifereetiiblee differential equation solvers are available nowadays, e.g. in diffrax. Finally, there is lot of interest in developing bioplausible learning algorithms where weights can be updated locally and in online fashion. Notable advances in this direction include [7]. To the best of our knowledge, none of these works cover the case of stochastic SNNs where the neuronal dynamics modeled as SDEs instead of ODEs.
58 46 40 This approach essentially analogous to duscretisetthen-pplee
3 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 SLIF neurons. Here the randomness appears both in the inter-spike dynamics well in the firing mechanism. To motivate the general definitions and concepts we start with an informal introduction of SSNNs.
3.1 Stochastic spiking neural networks
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 [3 and stochastic firing (50) 27 That is, the potential is modelled by eq. (T) Instead of firing exactly when the membrane potential hits set threshold, we model the spike times (event times) by an inhomogenous Poisson process with intensity X : T -> IR,, which is assumed
to be bounded by some constant C O. This can be phrased as an Event SDE (note that this is essentially the reparameterisation trick) by introducing the additional state variable s, satisfying
ds; A(vi-)dt, so log u
3

---

where u ~ On, 1). The neuron spikes whenever hits 0 from below at which point the membrane potential is reset to a resting level and we sample new initial condition for s, We can denote this first spike time by and repeat the procedure to generate sequence of spike times < < In practice, we reinitialize 8t at log u a for some a > O. It can then be shown that
(t < (x. min (1,exp a (a-
for t. E Tn+1).
It follows that Tn+1 Tn 2 a/C a.s., i.e. a 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 2 1 denote the total number of neurons. We model neuron & E [K] be the three dimensional vector y = (u*, ah )² the dynamic of which in between spikes is given by
du Mn (t ) dt to,B diç = --2 dt + oadB ds = > E)dl,
(2)
where B is a standard two-dimensional Brownian motion, o = (2), E R2x2, μ (11,112) E RR and E(= I- + R is an intensity function. As before, neuron & fires (or spikes) whenever s 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 will increment the input current of neuron by Wkj- Here w E R x is 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 O.
As presented here, there is way 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.
3.2 Model definition
Definition 3.1 (Event SDE), Let N E be the number of events. Let yo E IR be an initial condition.
Let p: I + Re and R + Rexd be the drift and diffusion vector fields. Let & IR - R and T : + R. be an event and transition function respectively. We say that (y, (Tw) is solution to the Event SDE parameterised by (yo, a. T,,,N) if yr u
N
1/2
Σ
De 1 (1)
On inf {1 --- E(11 0},
(3)
with E(y".) + 0 and
dyf = u(vi dt + ov) started at v yo: dyj u(v² dt + a (v." ddB, started at v". T (vr,
(4)
(5)
where B, is a d-dimensional Brownian motion and
(5) are Stratonovich SDEs.
In words, we initialize the system at yo. evolve using until the first time Tj at which an event
happens E1/2 ) O. We then transition the system according to y T (3 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
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 Event SDEs driven by continuous semimartingales has been studied in [31] Theorem 5.2] and (32). Under sufficient regularity of and a,, a unique solution to (4) exists.
We need the following additional assumptions:
Assumption 3.1. There exists e > 0 such that for all s E (0, 7) and a E im T it holds that
inf > s E(v) 0} > where ye is the solution to 4 started at y, a

---

Assumption 3.2. It holds that T(ker E) n 0.
Assumptions 3.1 and 3.2 simply ensure that an event cannot be triggered immediately upon tran- sitioning. This holds in most settings of interest. For example, for the usual deterministic LIF neuron im T 0 and ker & = 1 and the duration of the refractory period is directly linked to in Assumption 3.1
Theorem 3.1 (Theorem 5.2, 31). Under Assumptions 3.1
3.2 and with μ E Lip' and a E Lip" for 7 2, there exists unique solution (y, (Tn) 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, cases where the driving noise exhibits jumps. In the latter case, 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 unified treatment of differential equations driven by noise signals of arbitrary (ir- regulaaglr,,,,, and makes all proof's 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).
3.3 Backpropagation
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, p. a, and the event and transition functions &, TT might depend on some learnable parameters 8. We can move the parameters 0 of the Event RDE inside the initial condition you by augmenting the dynamics with the additional state variable 0r satisfying d0;- = 0 and 6,- = 8. Thus, as long as we can compute gradients with respect to yo. these will include gradients with respect to such parameters. We then require the gradients 0g J., if they exist. For this, we need to be able to compute the Jacobians dy = 0g Ve of the inter-event flows associated to the dynamics of y/f' and the derivatives Ör, == 0yo Tu. We assume that the event and transition functions & 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())) (T())a(y) = 0 for all E ker E Assumption 3.4. (v) 0 for all y ker E.
Assumption 3.5. VE(y)p(y) # for all y ker 3
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 holds for example if the event functions depend only on smooth part of the system. Assumption 3.3 what allows us to differentiate through the event transitions.
Theorem 3.2. Let Assumplions 3.1 3.5 be satisfied and (y, (Tn) the solution to the Event SDE parameterized by (yo, pt E,E, T. N). Then, almost surely, for any n E [N), the derivatives Sra and the Jacobians Oyf exist and admit the following recursive expressions
VE(UH. θT, VE(V"- ()((()))))
(6)
Oyf = (0pp, V/¹) 10 -VT -VT""-
θ.a].
(7)
where dyL and ÔTn are the total derivatives of sh. and Tn with respect to the initial condition yo- Oye yi denotes the partial derivative of the flow map of eq. (5) with respect to its initial condition, and VT E TRex and VE Lxe are the Jacobians 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) 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
5

---

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 (v( S and corresponding transition functions given by
Tk(v)= =
Th (v²))
where (z((²) = (v0, + Wkj, s²) if) j * le and (x") = (o+ Uresct: 1 log u - a) where Preset > i constant determining by what amount the membrane potential is reset. The addition of the constant a > 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 of exactly to 0 upon spiking would break Assumption 3.3 If one is interested in such a resetting mechanism it suffices to pick diffusion term 0 (v*) 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
3.4 Numerical solvers
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 0Øp, yi. We shall write SDE:Solveyyo- n1, a, 8, for generic choice of such a solver. Furthermore, if RootFind (yo f) is a differentiable root finding algorithm (here f (y,t) + R should be dif- ferentiable in both arguments and RootFind (yo. f) returns to E IR such that f(y. 0), then we can define a differentiable map E' yo + y by
RootFind (yo. ( o,
y" SDESolve (yo, μ. s,t,)).
Consequently, ExentSDESolve
a, E, N, can be implemented as subsequent com-
positions of T o E (see Algorithm [ This is a sescretise-then-opimmee approach (28).
2. Alternatively, one can use the formulas (6) and directly as 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 0ne, y/" with the derivatives of the given numerical solver. This approach is a solution between discretisethhen-ppiiiis and
3. Finally, one could apply the adjoint method (or optimise-ther ceeenise-teen--eee then-discretise as done for de- terministic 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-eee oottioot 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 B.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 time of the driving Brownian motion to 0 and picking step size sufficiently small.
Remark 3.4. In the context of SNNs, Algorithm iss 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 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 full SDESolve and additional RootFind each time an event occurs. This problem becomes especially prevalent for SSNNs with large number of neurons since in this case an event is triggered every time 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 2 0 and take this as the spike time. For the backward pass could then solve the adjoint equations (for which you need need to store the spike times 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 forward (or backward) pass of the entire SSNN scales as O(KS) where the cost of the forward (or backward) pass of single neuron and K is the number of neurons.
Algorithm EventSDESolve
N, Input yo: a, T, N, to, At,
1: y V/u
2: n + 0
3: + E(v)
4: while n < N and lo < < Tdo do
5: while
do
We assume for simplicity that
6:
yo x S.ESstvveStepww...oottt...
7:
8:
to + to + At
9:
e e E(y)
> Update value of event function
10:
end while
+ RootFind (No. E(SEESolvSSee (,, At)
> Find exact event time
12:
Un+1 < SSDESolvSStep a, Ae, A,,1+1)
> Compute state at event time & Apply transition function
+- T(v++1)
14:
+1+1
end while
Return (t))ngN,
4 Training stochastic spiking neural networks
4.1 A loss function based on signature kernels for càdlåg paths
To train SSNNs we will adopt 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 the introduction, classical signature kernels are not directly applicable to the setting of & 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.
The main idea goes as follows, a is càdlàg path, one can define the Mareus signature S(z) in the spirit of Marcus SDEs 42 43 as the signature of the Marcus interpolation of x. The general construction given in Appendix A The Marcus signature kernel is defined as the inner product k(r,, y) = (S(x), (()) of Marcus signatures S(x), Sly) of two càdlàg paths y. As stated in the first part of Theorem C.T 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)
v), E Ep..''ppk/2. (() x') 2E x²) Ey.s'' (y.3')))))
satisfies the property di (p, v)2 > μ for any compactly supported measures μ4,
Nonetheless, characteristicness ceases to hold when one considers measures on càdlàg paths that not compactly supported. In (8) the authors address this issue for continuous paths by using the
7

---

so-called robust signature. They introduce tensor normalization A ensuring that the range of the robust signature S remains bounded. The robust signature kernel is then defined the inner product ka (x. y) = (A S(x),, A o S(v))- This normalization can be applied analogously to the Marcus signature resulting in a robust Marcus signature kernel. In the second part of Theorem C we prove characteristicness of A for possibly non-compactly supported Borel measures on cådlàg paths. The resulting MMD is denoted by de.-
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 path- dependent wave equation [52]. Our experiments are compatible with both of these methods.
k==1 - sampled from underlying unknown target measure y we can train an Event SDE by matching the generated càdlåg trajec- tories {y') In ~ n using unbiased empirical estimator of dg (or dk, ) i.e. minimising over the parameters 0 of the Event SDE the following loss function
Given a collection of observed cådlág trajectories {>')
L
1
E jæi
k(x', .²)
2 mn
Σ
k(x',
1
Σ (M..v') 5
m(m 1)
n(n 1)
In the context of SSNNs, the observed and generated trajectories xt's and y correspond to spike trains, which are càdlàg paths of bounded variation.
4.2 Input current estimation
The first example is the simple problem of estimating the constant input current c >> 0 based on sample of spike trains in the single SLIF neuron model,
dup- ulc- we)dt + ad,,, ds. Xvr)d,
where A(v) = exp(5(u - 1), 1 15 and a varies. Throughout we fix the true c = 1.5 and set Preset 1.4 and 0.03. We run stochastic gradient descent for 1500 steps for two choices of the diffusion constant a. The loss function is the signature kernel MMD between a simulated batch and The test loss the mean absolute error between the first three average spike times. Results are given in Fig. For additional details regarding the experiments, we refer to Appendix[E]
the sample of spike trains |
In all cases backpropagation through Algorithm i able to learn the underlying input current after around 600 steps up small estimation error. In particular, the convergence is fastest for the largest sample size and the true e recovered for both levels of noise.
4.3 Synaptic weight estimation
Next we consider the problem of estimating the weight matrix in feed-forward SSNN with input dimension 4, 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 batch size of 128 and for sample size of 256, 512, and 1024 respectively. Learning rate is decreased from 0,003 to 0.001 after 1000 steps. The results able to reach a test loss of practically O. that is, samples from the learned model and the underlying model are more or less indistinguishabl.. 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.
given in Fig. 2 in Appendix B For a sample size of 512 and 1024 we
4.4 Online learning
In the case of SSNNs, equations ( Iead to a formula for the forward sensitivity where any propagation of gradients between neurons only happens spike times and only between connected neurons (see Proposition D.I.. Since the forward sensitivities are computed forward time together with the solution of the SNN, gradients can be updated online new data appears. As a result, between spikes of pre-synaptic neurons, we can update the gradient flow of the membrane potential
FFor simplicity we only compute approximation of the true MMD by truncating the signatures at depth 3 and taking the average across the batch/sample size.

---

0.25
0.25
Sumple sizo
32
1.2
1.0
0.8
0.6
250
500
750
1000
1250
1500
0
250
500
750
1000
1230
1500
Sample sioe
61
1.2
0.5
0.6
2000
5000
750
1000
1250
1500
250
500
750
1000
1250
1500
Step
Step
Figure 1: Test loss and estimate across four sample sizes and for two levels of noise a. On the left: MAE for the three first average spike times a hold out test set. On the right: estimated value of at the current step.
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 K2 gradient flows (one for each weight in the network).
On the other hand, the adjacency matrix the weight matrix forms a directed acyclic graph (DAG), three-factor Hebbian learning rules like those in [57 (2 easily derived from Proposition D.I For simplicity, consider the SNN consisting of deterministic LIF neurons and let No denote the spike train of neuron k, i.... N is càdlåg path equal to the number of spikes of neuron & at time t. We let 7 (t) (or 7 for short) denote the last spike of neuron & before time t. We shall assume that the instantaneous loss function L, depends only on the most recent spike times 7 "...,7K. Then,
Sw, Li =0,
N1 M1(1
where at is the eligibility trace and the first term be viewed as a global modulator, that is, top-down learning signal propagating the error from the output neurons The eligibility trace satisfies
da!* μr ) dt
afth
Ureseta (1 m
dNd. dbf - + +dNY
where the d.N terms are to understood the Riemann- -temlannStieltjjes 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, particular, Chapter 5). For general network structures one can use the eligibility traces as proxies for the the true derivatives o 7
5 Conclusion
We introduced a mathematical framework based on rough path theory to model SSNNs SDEs exhibiting event discontinuities and driven by càdlàg rough paths. After identifying sufficient
'Note that in the case of stochastic SNNs this term is not necessarily well-defined since semi mantingatss are general not differentiable wrt. time.

---

conditions for differentiability of solution trajectories and event times, we obtained a recursive relation for the pathwise gradients 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 define general-purposs loss function on the space of cådlåg rough paths to train SSNNs where noise is present both the spike timing and the network's 's dynamics. Based on these results, we also provided an end-to-end autodifferentiable solver for SDEs with event discontinuities (Algorithm I 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 beyond the scope of the paper. Examples that require much higher number of neurons than the two examples already discussed will be hard handle with the peccetize-then-opiiiiz approach for the reasons given in Remark 3.5
We think there are still many interesting research directions left 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 differs from the usual approach with surrogate gradients even in the deterministic setting, questions remain on how these methods compare for training SNNs. Furthermore, would be interesting understand 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, remains to be seen whether the discussion in Section c.4] lead to a bio-plausible learning algorithm with comparable performance to state-of-the-ar backpropagation methods and implementable on neuromorphic hardware.
Acknowledgements. Christian Holberg gratefully acknowledges financial support from Novo Nordisk Foundation through Grant NNF200C0062958 and from Independent Research Fund Denmark Natural Sciences through Grant 9040-00215B.
References
[1] Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International conference on machine learning, pages 214-223. PMLR, 2017.
[2] Guillaune Bellec, Franz Scherr, Anand Subramoney, Elias Hajek, Darjan Salaj, Robert Leg- enstein, and Wolfgang Maass. A solution the learning dilemma for recurrent networks of spiking neurons. Nature communications, 11(1):3625, 2020.
[3] Alain Berlinet and Christine Thomas-Agnan. Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media, 2011.
[4] Thomas Cass and Cristopher Salvi. Lecture notes on rough paths and applications to machine learning. arXiv preprint arXiv:2404.06588 2024.
[5] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[6] Ricky TQ Chen, Brandon Amos, and Maximilian Nickel. Learning neural event functions for ordinary differential equations. In International Conference on Learning Representations, 2020.
[7] Ilya Chevyrev and Peter K Friz. Canonical rdes and general semimartingales as rough paths. The Annuls of Probability, 4701)4200-463 2019.
[8] Ily Chevyrev and Harald Oberhauser. Signature moments to characterize laws of stochastic processes. Journal of Machine Learning Research, 23(176): 1-42, 2022.
[9] Nicola Muca Cirone, Maud Lemercier, and Cristopher Salvi. Neural signature kernels as infinite- width-depth-limits of controlled resnets. In International Conference on Machine Learning, pages 2552555.,,, PMLR, 2023.
10

---

[10] Thomas Cochrane, Peter Foster, Varun Chhabra, Maud Lemercier, Terry Lyons, and Cristopher Salvi. Sk- tretree a systematic malware detection algorithm on streaming trees via the signature kernel. In 202/ IEEE international conference cyber security and resilience (CSR), pages 35-40. IEEE, 2021.
[11] Sebastien Corner, Corina Sandu, and Adrian Sandu. Modeling and sensitivity analysis method- ology for hybrid dynamical system. Nonlinear Analysis: Hybrid Systems, 31:19-40, 2019.
[12] Sebastien Corner, Adrian Sandu, and Corina Sandu. Adjoint sensitivity analysis of hybrid multibody dynamical systems. Multibody System Dynamics, 49:395-420, 2020.
[13] Christa Cuchiero, Francesca Primavera, and Sara Svaluto-Ferro. Universal approximation theorems for continuous functions of c' c\(ddl) aag paths and IX evy-type signature models, arXiv preprint arXiv:2208.02229, 2022.
[14] Adeline Fermanian, Terry Lyons, James Morrill, and Cristopher Salvi. New directions in the applications of rough path theory. IEEE BITS the Information Theory Magazine, 2023.
[15] Peter Friz and Atul Shekhar. General rough integration, lévy rough paths and lévy-kintchine- type formula. Annals of probability: An official journal of the Institute of Mathematical Statistics, 45(4):2707-2765, 2017.
[16] Peter K Friz and Martin Hairer. A course on rough paths. Springer, 2020.
[17] Peter K Friz and Nicolas Victoir. Multidimensional stochastic processes as rough paths: theory and applications, volume 120. Cambridge University Press, 2010.
[18] Peter K Friz and Huilin Zhang. Differential equations driven by rough paths with jumps. Journal of Differential Equations, 264(10)::226-630)) 2018
[19] Wulfram Gerstner and Werner Kistler. Spiking neuron models: Single neurons, populations, plasticity. Cambridge university press, 2002.
[20] Julia Gygax and Friedemann Zenke. Elucidating the theoretical underpinnings of surrogate gradient learning in spiking neural networks. arXiv preprint arXiv:2404.1499 14964, 2024.
[21] Thomas A Henzinger. The theory of hybrid automata. In Proceedings 11th Annual IEEE Symposium Logic in Computer Science, pages 278 278-292. IEEE, 1996.
[22] Melker Höglund, Emilio Ferrucci, Camilo Hernández, Aitor Muguruza Gonzalez, Cristopher Salvi, Leandro Beebbezzrtttt,,urt,, and Yufei Zhang. A neural rde approach for continuous- time mon-markovian stochastic control problems. In ICML Workshop on New Frontiers in Learning, Control, and Dynamical Systems, 2023.
[23] Zacharia Issa, Blanka Horvath, Maud Lemercier, and Cristopher Salvi. Non-adversarial training of neural sdes with signature kernel scores. Advances in Neural Information Processing Systems, 2023.
[24] Hyeryung Jang and Osvaldo Simeone. Multisample online learning for probabilistic spiking neural networks. IEEE Transactions on Neural Networks and Learning Systems, 33(52:2034- 2044, 2022.
[25] Junteng Jia and Austin R Benson. Neural jump stochastic differential equations. Advances in Neural Information Processing Systems, 32, 2019.
[26] Danilo Jimenez Rezende and Wulfram Gerstner. Stochastic variational learning in recurrent spiking networks. Frontiers in computational neuroscience, 8:38, 2014.
[27] Hiroshi Kajino. A differentiable point process with its application to spiking neural networks. In International Conference on Machine Learning, pages 5226 5235. PMLR, 2021.
[28] Patrick Kidger. On Neural Differential Equations. PhD thesis, University of Oxford, 2021.
[29] Patrick Kidger, James Foster, Xuechen Li, Harald Oberhauser, and Terry Lyons. Neural sdes as dnfiiit-dimensionn gans, arXiv preprint 4rXiv:2102.03657, 2021.

---

[30] Franz J Király and Harald Oberhauser. Kernels for sequentially ordered data. Journal of Machine Learning Research, 20(31):1-45, 2019.
[31] Jaroslav Krystul and HAP Blom. Generalised stochastic hybrid processes as strong solutions of stochastic differential equations. Hybridge report D, 2. 2005.
[32] Jaroslav Krystul, Henk AP Blom, and Arunabha Bagchi. Stochastic differential equations on hybrid state spaces. Stochastic Hybrid Systems, 24(15-45):170, 2006.
[33] Petr Lansky and Susanne Ditlevsen. A review of the methods for signal estimation in stochastic diffusion leaky integrate-ande aan-utte-an--iree neuronal models. Biological cybernetics, 99(4-5):2532663 2008.
[34] Jane H Lee, Saeid Haghighatshoar, and Amin Karbasi. Exact gradient computation for spiking neural networks via forward propagation. In International Conference on Artificial Intelligence and Statistics, pages 1812-1831. PMLR, 2023.
[35] Maud Lemercier, Cristopher Salvi, Theodoros Damoulas, Edwin Bonilla, and Terry Lyons. Distribution regression for sequential data. In International Conference on Artificial Intelligence and Statistics, pages 3754-3762. PMLR, 2021.
[36] Xuechen Li, Ting- Kam Leonard Wong, Ricky TQ Chen, and David K Duvenaud. Scalable gradients and variational inference for stochastic differential equations. In Symposium on Advances Approximate Bayesian Inference, pages 1-28. PMLR, 2020.
[37] John Lygeros and Maria Prandini. Stochastic hybrid systems: a powerful framework for complex, large scale applications. European Journal of Control, 16(6)5533554, 2010.
[38] Terry J Lyons. Differential equations driven by rough signals. Revista Matemática Iberoameri- cana, 14(2):215-310, 1998.
[39] Terry J Lyons, Michael Caruana, and Thierry Lévy. Differential equations driven by rough paths. Springer, 2007.
[40] Gehua Ma, Rui Yan, and Huajin Tang. Exploiting noise resource for computation and learning in spiking neural networks. Putterns, 4(10), 2023.
[41] Georg Manten, Cecilia Casolo, Emilio Ferrucci, Søren Wengel Mogensen, Cristopher Salvi, and Niki Kilbertus. Signature kernel conditional independence tests in causal discovery for stochastic processes. arXiv preprint arXiv:2402. 18477, 2024.
[42] Steven Mareus. Modeling and analysis of stochastic differential equations driven by point processes, IEEE Transactions Information theory, 24(2):164-172 1978.
[43] Steven Marcus. Modeling and approximation of stochastic differential equations driven by semimartingales. Stochastics: An International Journal of Probability and Stochastic Processes, 4(3):223-22 1981.
[44] James Morrill, Cristopher Salvi, Patrick Kidger, and James Foster. Neural rough differential equations for long time series. In International Conference on Machine Learning, pages 7829-7838. PMLR, 2021.
[45] Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, Bernhard Schölkopf, et al. Kernel mean embedding of distributions: A review and beyond. Foundations and Trends@ in Machine Learning, 10(1-2):1-141, 2017.
[46] Emre O Neftci, Hesham Mostafa, and Friedemann Zenke. Surrogate gradient learning in spiking neural networks:
Bringing the power of gradient- bbaed optimization to spiking neural networks. IEEE Signal Processing Magazine, 36(6):51 2019.
[47] Ali Pakniyat and Peter E Caines. On the stochastic minimum principle for hybrid systems. In 2016 IEEE 55th Conference on Decision and Control (CDC), pages 1139-1144. IEEE, 2016.
[48] Alexandre Pannier and Cristopher Salvi. A path-dependent pde solver based on signature kemeess. arXiv preprint arXiv:2403.11338 2024.
12

---

[49] Christian Pehle. Adjoint equations of spiking neural networks. PhD thesis, Universität Heidel- berg, 2021.
[50] Jean-Pascal Pfister, Taro Toyoizumi, David Barber, and Wulfram Gerstner. Optimal spike- timing-dependent plasticity for precise action potential firing in supervised learning. Neural computation, 18(6):1318- 1348,, 2006.
[51] Cristopher Salvi. Rough paths, kernels, differential equations and an algebra of functions on streams. PhD thesis, University of Oxford, 2021.
[52] Cristopher Salvi, Thomas Cass, James Foster, Terry Lyons, and W. Y. The signature kernel is the solution of Goursat PDE. SIAM Journal on Mathematics of Data Science, 3(3):873-899, 2021.
[53] Cristopher Salvi, Maud Lemercier, Thomas Cass, Edwin V Bonilla, Theodoros Damoulas, and Terry J Lyons. Siggpde: Scaling sparse gaussian processes sequential data. In International Conference on Machine Learning, pages 6233 6242. PMLR, 2021.
[54] Cristopher Salvi, Maud Lemercier, Chong Liu, Blanka Horvath, Theodoros Damoulas, and Terry Lyons. Higher order kemeel mean embeddings to capture filtrations of stochastic processes. Advances in Neural Information Processing Systems, 34: 16635-16647, 2021.
[55] Bernhard Schölkopf and Alexander Smola. Learning with kernels: support vector machines, regularization, opfimization, and beyond. MIT press, 2002.
[56] Timo Wunderlich and Christian Pehle. Event-based backpropagation can compute exact gradients for spiking neural networks. Scientific Reports, 11(1): 12829, 2021.
[57] Mingqing Xiao, Qingyan Meng, Zongpeng Zhang, Di He, and Zhouchen Lin. Online training through time for spiking neural networks. Advances in neural information processing systems, 35:20717-20730, 2022.
[58] Friedemann Zenke and Tim P Vogels. The remarkable robustness of surrogate gradient learning for instilling complex function in spiking neural networks. Neural computation, 33(4))899-925
2021.
13

---

Appendix
The appendix is structured as follows. Section 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 [.2] but in the setting of Event RDEs as well as some preliminary technical lemmas needed for the proof. Secvion[givvss a brief overview of the main concepts kernel learning and presents our results on Marcus signature kernels along with their proofs. Section [D] derives the forward sensitivities of SSNN. Finally. Section covvrrs all the technical details of the simulation experiments that were not discussed in the main body of the paper.
A Cådlåg rough paths
Marcus integration developed in 7 preserves the chain rule and thus serves as an analog to Stratonovich integration for semi semi-martingales with jump discontinuities. In particular, it allows to define 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 in 77.
Let C([O, TJ. E) and DI, E) be the space of continuous and càdlåg paths respectively on [0, with values in a metric space (E, d). For p 2 1, let (p,, T], E) and [1, TT, E) be the corresponding subspaces of paths with finite p-variation. For any N 2 Let Let G. (Rd) be the step-N free nilpotent Lie group over I2d endowed with the Carnot-Caratheodory metric d. Let n (R4) == ((0, 7T, Glpl (IR)) and D2 (R²)= == (10, T1, Girl (R²) be the space of weakly geometric continuous and càdlåg p-rough paths respectively with the homogeneous p-variation metric
yx = max sup 1skštpl DC
d( di
Define the log-linear path function
Φ G. (Rd) x G (R2) - C(10, 1], G.¹ G( (R4))
ba) + exp((1 - log + log b).
where log and exp are the (truncated) tensor logarithm and exponential maps on G. (Rd). If N = 1, then GO (R2) z IR + d and ó(a, b). = (1, (1 - t)a + tb) is a straight line connecting a to b in unit time. For any x E [ TJ, G. (R))) we can construct a continuous path x E C(10, T), GN (R')) by adding fictitious time and interpolating through the jumps using the log-linear path function according to the following definition.
(N, (R')., let 71, T2: Tm Definition A.1 (Marcus interpolation). Let Let 2 1. For D( T1, GN
be the jump times of x ordered such that ((X,,,,x,)) 2 d(x;,_,,,,) 2 where 0= < I do is the number of jumps. Let (rk) be sequence of positive scalars T > 0 such
2 d(x X,m) )
that r = Ek=1 Tk < +00. Define the discontinuous reparameterisation rj [0, T] - [0, T + r] by
((( = 1+
Th (mS())
k=1
The Marcus augmentation x E C((0, T + r), N. (R) of x is the path
X1 10(x Xyx (s (1-
iʃs = (1) (() some 1 E |0, T],
x =
if [nT(T=) n(Tx)) 1 S k < + 1.
The Marcus interpolation & C( N. (Rd)) )} x is the path x M Tr where 1y (t)
t(T +T/ is a reparameterisation from (0, T to [0,T T + r). We can recover x from & via x = x 7/x by considering the reparameteriisation 7/x n, o 17.
Once the Marcus interpolation is defined we can state what we mean by solution to differential equation driven by a geometric càdlàg rough path.
14

---

Definition A.2 (Marcus RDE), Letx E D (R'1) and = (fi..... fa be Lip" vector fields on R
with y > p. For an initial condition a E R", let y E (1, be the solution to the classical RDE driven by the Marcus interpolation & E 22 (IR))
dýt (ū))d
Define the solution y Dp([0, T). R to the Marcus RDE
dy: (3)) o dxi,
(8)
to be - o The where is the reparameterisation introduced in Definition A.1
A.1 Metries on the space of càdlåg rough paths
Chevyrev and Friz 7 introduce metric ap on 22 (RR) with respect to which 1) geometric cådlåg rough paths can be approximated with sequence of continuous paths [ Section 3.2] and 2) the solution map (Yo. x) (x, y) of the Mareus RDE (8) locally Lipschitz continuous 7 Theorem 3.13].
We write A for the set of increasing bijections from [) to itself. For À A let [|] suPre(o, [A(t) - t). We first define the Skorokhod metric as well as Skorokhod version of the usual pariaatio metric.
Definition A.3, For p 2 1 and x, E D2([0, we define
[A]. sup d((x J tEjU,Ty
xx (x, = inf max AEA
(x(x,y) = inf max (xo
AEλ
turns out that the topology induced by op 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
E 2. (R*) and f (fi...., fa) family of vector fields in (R²) with y p. let y s, x)) denote the solution to the Mareus RDE dy: f(y) dx, initialized ys y and evaluated at time t. We define the set
(a', U/) a, J, = ((, (a,,
0:::0,:1000000 1; -b}
and, on it, the path function
(,((a, (a', ()))))))))))))))))))()))))))(( (1, a():))
Finally, let ((0, 7), GlPl (Rd) R be the space of càdlåg paths (x, y) on (IR)(R x R of bounded p-variation such that (z-,zz, 2t) )) Jy for all jump times of keep notation simple, we shall write DI when this does not cause any confusion. Naturally, if y is the solution to the Marcus RDE dye = f(v) o dx,, we have (x,y) y) E DJ. For a z = (x, y) E D we may define the Marcus interpolation by interpolating the jumps using Φj. Let 2 denote this interpolation but with Te replaced by or% for 8 > 0 and similarly for x with x E 2227 (Rd)).
Definition A.4. For f (f....., fa) a family of vector fields in z, E D with (x,9) and x' (x's) and define
with y > let
a,(x, x') = lim op rx( x'd), 8
0 12 z') = lim op pp 225).
6-+0
Remark A.1. It is proven in 77] that in both cases the limit in ap exists, is independent of the choice of rk. and that it is indeed a metric on 22 (IRd) resp. D/.
Theorem A.1 (Theorem 3.13 Proposition 3.18, 77. Let S (fr,..,, .S)) be a family of vector fields in Lip" Lip"-((RR") with 7 > p. Then,
15

---

1. The solution map
R (02(R")) ap) - (Dp. ap) (yo, x) 7 = (x,
of the Marcus RDE dy.. [wv) o dx, initialized at E Ik is locally Lipschitz.
2. On sets of bounded p-variation, the solution map
R x (22p (R2)) 0 (D,(10, 1),
(yo. (6.x))+9 y
of the Marcus RDE dy: f(ye o dx, initialized at yo E R2 continuous.
Now, let Cj (IR)) be the space of absolutely continuous functions on R.. Definition A.5. We define the space of geometric càdlàg p-rough paths 22 (R²) as the closure of Ci. (RR) in 22 (IR") under the metric 0xp-
Remark A.2. A càdlàg semi-martingale a E (,( 7J, 12)) can be canonically lifted to geometric càdlàg p-rough path, with p E 3), by enhancing it with its two-fold iterated Marcus integrals, i.e.
x 25,t1 /
F.)) @ odx.) (IR))
where the integral is defined in a similar spirit to Definition A.2 (see, for example, [8] 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 æ has a.s. continuous sample paths, agrees a.s. with the solution to the Stratonovich SDE. See, e.g., Proposition 4.16 in 77.
A.2 Signature
The extended tensor algebra over Rd is given by
T ( (TR*)
= II
(R= 8n
n=0
equipped with the usual addition + and tensor multiplication @ An element T ((R)) is formal series of tensors
(a", a .!....) such that a" E (Re))or. We define the projections Tn T ((Id) + (I2))n given by n, (a) = a". Let T((T)) be the subset of T((R*) such that the (a) = 1 for all a E (((R))) Finally, we define the set of group-like elements,
((+) = {a a I ((x4) (a(a)) E N.N (IR") for all n 2 0}
Definition A.6. Let p and E 2.P (R.). The signature of is the path S(x) [0,7] 1 G such that, for each N 2 0,
dS(x), = S(x); @ odx,, S(x) = 1 G. (Rd).
(9)
Remark A.3. Uniqueness and existence of the signature follow from the continuous analog. Indeed, by definition, (9) is equivalentt to a continuous linear RDE.
Remark A.4. The signature, as defined here, is also known the minimal jump extension of x and was first introduced in [15]. It was further explored in [3 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 to tree-like equiv- alence. Two continuous paths y,y are said to be tree-like equivalent if there exists continuous non-negative map h: [0, T + R4, such that h()) hTT and
X,, Ys.ell & h(s) + h(t)) - 2 inf h(u). aE|s.t)
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, are tree-like equivalent. It is straightforward to check that this indeed an equivalence relation on (Rd). Perhaps more interestingly, we obtain the following result. For ease of notation we shall henceforth mean S(x)r when omitting the subscript from the signature.
~i y. if their their Marcus interpolations (see Def. A.I). & and y.
16

---

Proposition A.1. Letp 2 1. The map S(:) 2 (B<)) - G is injective up to tree-like equivalence, i.e, S(x) = S(y) iff x -xxxx y.
Proof.
The result follows from the continuous case upon realizing that S(x) 5(x) and analogously yory.
A.3 Young pairing
In many cases, given a geometric càdlåg rough path x (RR with p [2,3) and a path h E Di((O, Ti, RR) of bounded variation one is interested in constructing new rough path y ( (Rd+t) such that the first level of y is given by y (z, h). In the continuous case this can be done by using the level two information x² and dh dh to fill in the corresponding terms in y² 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 h2 = t. that is, we want to augment the rough path with an added time coordinate (see Def. A In the cådlág case one needs be more careful in defining the appropriate Marcus lift.
Definition A.7 (Definition 3.21 Z. Letx E D (R") (B*) with p 2 1 and E D ([O, R*). Define the path (x, h) and the corresponding Marcus lift 2 (2, The Young pairing of x and is the p-rough path h(x,h) E Sp (Ra-c) such that
h) P(&,h) Tl
where P(&, hi is the usual Young pairing of 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 he = t. It turns out that this pairing is continuous as a map from (RR.).((.)) ) to S..p
Definition A.8. Ler x DP (R*). The time augmented version of x is the unique rough path & D. (R(+1) obtained by the Young pairing P(x, h) of with the continuous bounded variation path h4 = f.
Proposition A.2. Let E (1,3). Then, the map x +> & is continuous and injective as a map from Do (RR) Dp (Ra+1).
Proof. Let === = PP,, be a metric space when equipped with Og-- Fix x E X' and let xt be a sequence of absolutely continuous paths converging in d' to We shall first show that H" then converges to x. Since x" does not have any jumps and any reparameterisation of x. is still absolutely continuous, we may assume that
ax(xx
lim dy 8
dx (x* x²) +
for n co. Define (x, h) and add = =x5 1 the Mareus interpolation with 1x,5 the reparam- eterisation such that 2 1x.55 Furthermore, let P(x, h) be the Young pairing of and h. By definition,
0g (() .N) lim op 8-+0
(((* RR)) P(a", h))
S lim dp 8-1
(P(*' h:²), P(x", h))
lim C 6-+
(d,((< "") di di(ht (h. .)) +
for n where C is just some generic constant depending only on p. The last inequality follows
from 9.32 in (17). Thus, y E X' is such that y)) < e, we can choose another sequence y" of absolutely continuous paths and N 2 large enough so that
ap (P(x, b), Ply, h)) < 2e ap(P(a" (h), P(y", P(y",h)).
17

---

aplx".'' S 2e
for all n 2 N. By Remark 3.6 in 77, we then have that, up to choosing a large N, dz(x", (x",,") y²a) S for all n 2 N and therefore, once more appealing to Theorem 9.32 in [17]
Op (P(x", h), P(y", h)) S 2Ce.
In conclusion, ap (P(x, P(y, - 2(1 + C')e which proves the result.
Injectivity follows from [13..
A.4 Event RDEs
The results of Section 3.3 hold in more generality. In fact, we define Event RDEs similar to Definition 3.1 where the inter-event dynamics are given by Marcus RDEs driven by càdlag rough paths. Utilizing the correspondence between solutions Mareus 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 2 and N
be the number of events. Let x E D (R") and f (fi...., fa) be a family of Lip" on Iº with y > p. Let E R + and T E. + R be an event and transition function respectively. We say that (y. (m)))i)) is solution to the Event RDE parameterised by (yo, x, S.E. if vr.
N
V/ Σ
yi 1 (mitntt)
Tn inf {t Tn-1 -(v(--) = 0},
(10)
a=0
with E(v,)) ) + 0 and
dy: = f(()) dxt. dy; (w²) dx started at v"r. T (V---).
started at 30 = you
(11)
(12)
Existence 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 exists. In fact, the solution map ys x (s, t) + ye is diffeomorphism for every fixed (see, e.g., Theorem 3.13 in 77).. It follows that we can iteratively define a unique
0 S S < sequence of solutions y" Dp([t,, 77, R"). Finally, rough path x is the Marcus lift of a semi-mann martingale, the inter-event solutions agree almost surely with the solutions the corresponding Marcus SDE.
mentioned in Remark A.2 if the driving
Theorem A.2. Under Assumptions 3.1 3.2 there exists unique solution (y, (Tn) w to the Event is the Marcus lift of a Brownian motion, the solution
RDE of Definition A.9 Furthermore, coincides almost surely with the solution to the corresponding Event SDE as given in Def. 3..
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 much larger family of driving noises such as fractional Brownian motion or even smooth controls. Also, since the driving rough path allowed to be càdlàg, the model class given by Def. A.9 includes cases where the inter-event dynamics are given by Mareus SDEs driven by general semi-martingales.
B Proof of Theorem 3.2
The proof Theorem 3.2 presented below covers the case where (y, (Tn) is the solution an Event RDE. Throughout we consider vector fields μ E Lip", a Lip"+e and specialise Event RDEs where the inter-event dynamics are given by
dyj u(y²" )tt + (vy") dx4,
(13)
where x E (PP(() The notation above deserves some clarification. One can define the vector field f= (1u, .) and the Young pairing x, of and ht l. Assuming μ E Lip2+e then view y/ as the unique solution to the Marcus RDE
dy f(v:) x
18

---

Alternatively, if one is not ready to impose the added regularity on the drift μa one can view [3] as a RDE with drift in Ch. 12 in [T7]. 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 semi-martingale, we shall need the following two additional assumptions:
Assumption B.1. For any n E [N], there exists non-empty interval I, = (Tn on In + 8x) such that x is continuous over I,. In other words, the càdlàg rough path x, does not jump in small intervals around the event times (Fn)-
Assumption B.2. For all 0 S S N we define Tn - -6x/2 and t, Tx+1 + 1/21/2. It holds that x E 00.p ([sn, f,), ). i.e., x is geometric p-rough path on the intervals [s,. tm] Remark B.I. Note that Assumption B. trivially holds if x is continuous. Otherwise, it is enough to assume, e.g., that x the Marcus lift of finite activity Lèvy process. Furthermore, by the properties of the metric ap, if x is the canonical Marcus lift of a semi- semi-martingale I E Dy([s, t),
[), Rd-1), then there exists sequence (x" of piece-wise linear paths xm E ([0, T}, 1dd such that
ap.(snat]] (xx) x)
as m 0C a.s.
See, c.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.
We shall need technical lemmas for the proof of 3.2
Lemma B.I. Assume that Assumptions 3.7 3.5 and
B.
B. B.2 are satisfied. Then, there exists open
bull Bo such that the following holds:
1. For all a E Bo, [T(a]] = N.
2. For any [N), the maps
Bo 3 Dam (T(())
are continuous.
3. For the sequence (x.") as given in Assumption B.2 and (y"T, (TT)) w the corresponding Event RDE solution, for all n E [N] it holds that
lim sup (17m" (a) - T, ra((a) + (T (a) (a)Ta(aa)) = O. 00 oBBo
Proof. Recall that (y, 8, 1; x) is the solution map or flow the differential equation
dy (yJ) o dx.., y/
y
evaluated at time The first step will be to prove continuity at Yo- In particular, let yū" E O approach yo for m going to infinity and denote the solutions the corresponding Event RDEs by yn (T > We claim that limm Nm = N and
lim m-
To Tnx
lim m00
man1
V
To this, note that, by Theorem A. there exists a sequence Am A of continuous reparameterisa- tions such that Aml + 0 and
sup (s.0)AT
(P s, s,);x) ( At" +
(14)
for m + oc. Note, furthermore, that (vv"" t; x Xm) = ( Am (s), Am (f): for all (s, t) e
Ar. We let fiTh ( (
No..
be the solution to the Event RDE where (yo, x) is replaced by (vi" x
Am). It suffices to prove that, for all < n N.
lim 4m
Tu
lim
m
Vr.
(15)
m
m00
Indeed, since pm = A-1 (T. and Am + 0. it then follows that ,. - Tn for m - c0. Furthermore,
we have 7Tm
-m.n- mnn
7m
19

---

We shall proof 15
using an inductive argument. We have that
y (yo. 0, x;
E [0,
y/ m.0 ( 0, 1 Am), Vt E [
Now fix some 0 <6 < d1 where ón is given in Assumption B.1
Note that [E(')) > 0 for all N such that, for all m mp.
t E (0,71- -(] and therefore, by [4), follows that there exists an mg
inf rE|0,T1-4
(x( 0, xx Xm)) 0
so that fi 2 -. Next, for some small 0 < 1) < e, Assumption 3.4 and the Mean Value Theorem imply the existence of oj = -r v (1 nn ) +op and an rn v + (1 raj )y with en Ta E (0, 1) such that
8 (w = = VE(a,)
ply))dye:
& (wh-w) = E (v%) -
But then, by Assumption 3.5 and the fact that E1,², = 0, for 1 small enough, ((r" and Ev -n) must lie on different sides of O. Assumption B.1 and eq. (14) then yield the existence of a m2 2 mo
such that 7.mh < + 7 S T and in +o] [0t m. 0 for all 2 mi. It follows that Ti + T1. Finally, note that
yo
+
Another application of (14) shows that the first term the right hand side goes to 0 for and second term vanishes by Assumption B.1
To prove the inductive step, assume that 15 holds for i S n. For all t E [Tn, In+1] it holds that mnn Vt
1b TT
-m.n mm V=
Tn tix (exx vi² + (T (V²-2) TT,,;
- and, since Fem
y 2-1 Tm -> Tps and T is continuous,
lim (E(T..T]
sup
① ( (Dmm.
A t: Xm) +b (T (T(?TT) Tny Tn, xx
whence the same argument as above proves that 15 also holds for n - 1. This completes the proof of the claim.
Now, by continuity yo, it follows that there exists some small r >> 0 such that for all a (y0) it holds that [raa = N and Tr (a) E (Tn 5n/2 Tn. + /,) for all n E [N] where 8, is as in Assumption B.1 Furthermore, since Assumption 3.1 3.5 and B. B.2 still hold for (yo), the
same argument as above can be applied to show that Ta (a) and y "r,(a) (a) are continuous at a. This
proves parts] and 2
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' not depending on æ such that
(p (a), 8d (a)) S Cap.jo.t) (x'",x) ,))
Since the latter term does not depend on y and goes to 0 for m going to infinity, we find that
lim sup (p.[0.ta] (m))) (a). m=+00 (20)
0.
(16)
Recall, (a) is the continuous path obtained by the Marcus interpolation with ora instead of rk and similarly for mm) Note that ym (a) = (a) by continuity. Letting - (a) and (a) denote the first event time of (a) and (a) respectively, we have, for all m E N
sup (a) TI(a)1 < ac B,
sup lim (1 - (a) + (+((a) - -())))
WE B.(o)) 8-0
20

---

for forall a E Bo and x is continuous over this interval, it follows that (f(a) - [=(a] goes to 0 as - 0 for each a E Ro- Furthermore, by definition of the metric Opp eq. (66), and the fact that y/
Now, let Bu = (y0)). Since [1) E (71 (
(a) = a = va(a), for each argument the one employed in the beginning of the proof then shows that 177" =f(a) - as
Bo, similar
5 > O and, thus, limm SuPaEBo [T] [a)a = 0. Finally, starting from the inequality
m.0 (a) -100--ane (a) - (a (I) -f((a) (a) + (T (r)
6.0
(a)
and taking the limit as -> O and then the supremum over x E Bo on both sides, we can argue in
exactly the same way to show that holds for 1. We can then argue by induction, just as in
the first part of the proof, to show that holds for all subsequent event times as well. Thus, the set Bo satisfies all the stated requirements.
O
Lemma B.2. Let Assumption B. hold and xmth be as in Assumption B.2
Then, for all n E [N] and
p' > p.
lim dp-1ee (x x) = 0, for any Tn 5n/2 2 S 8 < 10000
< L Tn + 8n
Proof. Fix some n E [N], p' > p and Tn - 55n/2 <<<<t>> Stt Tn + 8n Note that, for any continuous N. and 5 > 0, it holds that
reparameterisation À E A,
de/,e,,((""" ,x) S d'' d''.[s.0 (x" ,x A) (on,t,) X, x,x%) + 1e., ,x),
where x is the Marcus interpolation of x over the interval [s,, tnl. Taking the infimum over À A and the limit as 8 + 0 on both sides, we obtain
(s. (x" x) S Ctp', (sn (x. x) + lim dyy,[ssa (Se), x). 8-0
The first term on the right hand side goes to 0 as Assumption B.1 x is continuous on (T,, 6n Tn + 6,). it follows that [o.. (8.. x) goes to 0 for
by Assumption B.2 Furthermore, since, by
5 20. But the result then follows from Proposition 8.15 and Lemma 8. 8.16 in 17).
L
Proof of Theorem 3.2 Step 1: Assume that x E Ci C1 ((=)))) By 17 Theorem 4.4]. the
Jacobian Qyß exists and satisfies 7 for all t E [0. T.). We shall prove that relations (6) and (7) hold for all n E [N] by induction. Thus, assume that oy, and Ôr, exist for all / E [Tk, 7k++) and k S and satisfy the stated relations. To emphasise the dependence on the initial condition, we
will sometimes use the notation y" = (20) and Tn (n())) for the solution of the Event RDE
started at you We want to show that, for arbitrary E IR the following limits
yi n. ya
Tn and lim c
lim
for t E [Tn, Tn+1)
6
6
exist and satisfy the stated expressions, where 75 Ty (yo + he) and y y + he)-
For any E >> 0, because È is continuously diffcrentiable, the Mean Value Theorem implies that there
exists Ce E IR on the line connecting v²,,- to v"r and another de E R on the line connecting n-1,t to ش such that
& (yrh (0"-1) = n-
+ VE(c.) v²- - -V7a )
Urg
n-1
+ VÊ (c+) (K(V)- ((m - -ii) + a (V7-1) (or. - -xrz) + o([TT 1n1))
Uri (
+ VE(C.) 0 (v- - - 1111
=
+ VE(( (e (0w".!) ha o(t))
=
where the last equality follows from the induction hypothesis. We have E(vZ,1) = 0 = E(1," Thus, by rearranging, we find that
(0(W,-10097 1h
+ o(1)
VE(YR-1) (M(UA.)) a83. 1) ) 2m Bee To -
VE(V--1)08) + o(1)
VE(HH- ")((
21

---

where the second equality follows from Assumptions 3.4 and 3.5
Assume for that v n-- By another application of the Mean Value Theorem, there exists E IR on the line connecting ver!! to Vra such that
TT Win V T (yTT))
V"= T
U2: n-1 (V - VT
The + =M(=): (Tn - = (y=i1) (IT,, - -rrr)) -T V//
- -VT(ce) (u(V ((m - T + (V%1 (I+, -b+z) + +o((Tn r,)
= =T (v+E - -
U/r4 (u(()=;)) - VT(c, (u(()) (Tn - T,
+ (o(VC.) (V(c)) (ITn - -xr;) + ol(T - -1l)
Therefore Vhn
V
(u(v²,)) VT(" ( θr + o(1)
where used Assumption 3.3 the chain rule and the existence of θTn Finally, for any E (Tay Tnt1l equation (7) follows from the fact that we can write y? = ( 8, (. for all Tn S 8 < t. In particular, by the chain rule, find that
Oyf 0ay ((v s, s,ty0 S=TH Oze, yi [0y"l
Step 2: Consider now the general case of x (P") and let (3"", (T ) denote the solution to the Event RDE where x is replaced by the piece-wise linear approximation x With oy and 0Tm denoting the corresponding derivatives, we saw in the previous step that both exist and satisfy (7. We let R2 and Pn denote the right hand side of (7 and (6) respectively. This step consists of proving that, for n E [N] and E (Tn)tt).
lim {|T - m0-500
mal + 1Bte th
(17)
and for some open ball Bg around yo
lim sup { (10Ti - px((a) (c)- m- BB
+ 10yg "0y....000 R. (a)]l =.
(18)
By Lemma B.1 and continuity of T we have that Tyy m,n-11 converges to T(y"-1) and T converges to Ta as m + +o0. Then, because yi = (T((""" Tn, Tny 4; x) and ye m.n
(T (""..")), ), rm t; t;x"),, equation (17) follows from, Lemma B.2 and Corollary 11.16 in 177 In fact, since Bo was constructed in Lemma B.1 in such a way that Tn (a) tn for all a Bo we also get that
lim sup ==000 aBBn
(VF.,.))))( [n))), (.)) - i min (a) 1
(1 m.n (1 (a), Th (a).t;z")) =0.
by the same corollary in Thus, to prove 18) it suffices to show that, for all n E f),..,,N,
lim sup m 000 aBB
Oy To (a) - R (a) =
We shall prove it using another inductive argument starting with n 1. In this case it suffices to
show that
lim sup 10,4 (a, 0, (a): 0 4 (a, 0, -|" (a): 0. ac
By 77] Theorem 3.3] we know that the above holds if (a) and 7 (a) are replaced by t. + 61/2. Now let --! be the reverse of the flow map 4 that is,
(a1, s,t;)) as + b(a0, 8, t:x) 01.
22

---

From Lemma B.1 it follows that o (a) (a) = Φ fo:)) and, for m large enough,
lui x") But the result then follows from Lemma B.2 and 17 Corollary 11.16]. To prove the inductive step, assume that (8) holds for all i S n I. Again, by inspecting 6 and 7 and using the inductive assumption, one finds that it is enough to show that
m.0
(a) = + Ueo m.0 (T"
lim sup mh 720 BB
0 Ov 1, 2 (V 1n Twi - I
x=" =0,
where we suppressed the dependence on a for notational simplicity. This is done exactly for y and completes the proof of Step 2.
Step The third and final step is combine Step I and finish the proof. So far we have proven that 1) the theorem holds for continuous paths of bounded variation and 2) (r : ") converges to (Ta, VE') and (Or, 0y, m. (a)) converges uniformly to (Pn(a) R2 (a)) over a E Bo for all / E (Ta, (n) and E [N]. From these results it immediately follows that (Ta(a), Vi (a)) is differentiable at a = yo with derivatives given by (Pn (yo), Rg (y()) for all &. (Tn: (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 Tn <s < Tn + 6n/2 < t < Ta+1,
yi (v55 8, xx = Dy Ogr D(y s, f Rx (yo) R² (J6).
C Kernel methods
We give here a brief outline of some of the most central concepts related to kernel methods. For more in-depth introduction we refer the reader to [45 55 3 Let X' be topological space. We shall in this paper only be concerned with positive definite kernels, that is, symmetric functions k X' X + RR for which the Gram matrix is positive definite. To such a kernel one may associate feature map X' + I2 such that x + k2 k(zz,-), ). A reproducing kernel Hilbert space (RKHS) is a Hilbert space H R- such that the evaluation functionals, eV2- : f + fax, are bounded for each E X. For all positive definite kernels there is unique RKHS H. C R2 such that f(x) (kz, f) for all f H and x E X. This is also known as the reproducing property. Furthermore, with H denoting the linear span of {kz I E }], it holds that H = H. i.e., H is dense in H. Two important properties of kernels are characteristiceess and universality.
Definition C.I. Let k : X X + TR be a positive definite kernel. Denote by H the linear span of (k2 x E X} and let F C xR. be topological vector space containing H and such that the inclusion map i И + continuous.
* We say that k is universal to F if the embedding of L: H + F is dense.
We say that k is characteristic to F' the embedding μ. F + H', D + D/n 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' which case k being characteristic implies that the kernel mean embedding μ +
Ex= (-) is injective.
Remark C.2. Often times, instead of starting with the kernel function k and then obtaining the RKHS, one starts with feature map F X' + H into a RKHS and then defines the kernel as the inner product in that Hilbert space, i.e., k y) (F(x), Fly)ru. 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 Definition 6 in [8] for precise statement.
C.1 Marcus signature kernel
The definition of the signature kernel requires an initial algebraic setup. Let (-,-) be the Euclidean inner product on Rdd Denote by & the standard outer product of vector spaces. For any n E N. we denote by (:, >n, on (Rd)) the canonical Hilbert-Schmidt inner product defined for any a (a1...., .,) and b (b1...., .,,) in (Rdy@n as (a, b),, I[ (a,,bi), b-) The inner product (:
23

---

on (R*))" can then be extended by linearity to an inner product (:.-) on T((R)) defined for any (1, (1....... and b = (1, ..... in T((T)) as (a. b) = + Em=1 (an; bn)nn
To begin with, let X' = Dj (10, T7, Rd). If x E X' is càdlăg path, we can define the Marcus signature in the spirit of Mareus SDEs 42 4 as the signature of the Marcus interpolation of x. This interpolation, denoted by â. is the continuous path on [0, T' obtained from I by linearly traversing the jumps of z over added fictitious time r > 0 and then reparameterising so that the path runs over [0, 2"] instead of [0, T' + r]. The general construction is given in Appendix A If z is continuous, and å coincide; thus, without any ambiguity, we can define the Marcus signature (ax) of general bounded variation cådlåg path as the tensor series described above, but replacing x with è (see also the definition A.2
Since the signature is invariant to certain reparameterisations (Proposition A.I D. it is not an injective map. Injectivity is crucial property required ensure characteristicness of the resulting signature kernel that we will introduce next. One way of overcome this issue is to to augment a path æ with a time coordinate resulting in the path = (z, t D The Marcus signature kemel is then naturally defined as the map k : X x X - I such that k(zz, y) (S( S(ū) for any I, y E X. As stated in Theorem C.I] this kernel is universal on compact subsets K C X and, equivalently, characteristic to the space of regular Borel measures on K. However, these properties do not generalize to the whole space (xY, R) of bounded continuous functions from X' to R.
In [ the authors address this issue in the case of continuous paths by introducing the so-called robust signature. They define a tensor normalization continuous injective map
A: T ((R") {a T ((R))
||||| S R
for some R > 0 and such that A(a) = (a", X(a)a1, X(a)aa,....) for some A: 7 ((R4) - oo).
Now, let p E [1, 3) and take (Rd) to be the space of absolutely continuous functions on R Recall that 22 (IE") is the closure of Cj (IR)) in (R)d) under the metric ap. Throughout we let X' = -6, (IR5) be metric space equipped with ap. Naturally, we can then define the signature kernel on X' by k(x, y)
( (()) and, similarly, the robust signature kernel hA yx = (A))())), (y)())) where A is tensor normalisation.
Theorem C.1. Letp 2 L, A a tensor normalization, and K X compact under ap. Then,
(i) The signature kernel l is universal to F = C(K, TR) equipped with the uniform topology and characteristic to the dual F', the space of regular Borel measures on K.
(ii) The robust signature kernel KA is universal to F Ca(X. 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 [3]. For part (ii) we shall proof that the feature map F A o S is universal and characteristic. The result then follows from Proposition 29 in [8] We start by defining P = X/ ~ where the equivalence relation ~ is defined in Appendix A.2 We equip P with the topology induced by the embedding
S P + [((2+))). By Proposition A.1 F is continuous and injective map from P into bounded subset of (((Ra+))). Thus, 7 = ((E, F)) & T(()"+)))) is 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 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 de denoting the MMD for given kernel k X + R, the following is a direct consequence of Theorem C.I
Corollary C.I. Let 2 1, A a tensor normalization, and K C X compact under Og. Then, de is metric on M(K) and de, is a metric on M(X).
'If x is a càdlág rough path, this is done via Young pairing which results in a càdlàg p-rough path. &, where the first level is given by (æx, For more information, we refer Appendix A.3
24

---

D Forward sensitivities for SLIF network
In the general SSNN model, Theorem 3.2 gives the following result.
Proposition D.1. Fix some weight wij E w. a neuron k E [K] and let G denote the gradient of (o*, (²) wrt. Wij time t. Furthermore, define 7 {0, + IK² such that To [MI, ( -2) 71 (H1, Uresstt and let R22x2 be the drift matrix in the inter-spike SDE of ( Then.
g = (it=s) (gt You Ow., s dadgre2).
(19)
where en B2 is the n 'th unit vector, i is the neuron in Paye U [k} with the most recent spike time before t, and we denote this spike time by s. If is a spike time of neuron it therefore follows that
Mu
Jusy Eprev
Sigrer VA(ut)eTTT dr
Sw.tt
(20)
Maxt)
where Eyreu is the previous spike time of neuron k. In the case of a deterministic SNN, formula (20) is replaced by
Oun t == of
(21)
pi (2 x
Proof. Throughout fix some > 0 and let < denote most recent event time preceding with ( the index the neuron firing at time s. We define the process dw, = Odf with wo = Wij and with a slight abuse of notation we shall write yt (uî at, sf wi). We will leave out the event index a for notational simplicity. Since v depends on y,, only through v and VT, block diagonal, direct consequence of eq. 7 is
g = (I 0) Qy: (VT) (m.-))mmmmm vē_ (m(vt) VT (.::))(vt)))) Smy.))
where i,,,,,,, 8, w) (Mi - -2, A))), If E Pac U (k), then T id and therefore Of = (1 0) 0g Om,, V²- One can then reapply the formula above until 1 E Pac U [k]. By the flow property, it follows that we may assume without loss of generality that / E Paι U (k). This leaves us with two cases. We define 2 = (of +t) so that (7 0) dygut = d.42f and aw., 2 = G.. Furthermore, let a = 5z155k.
Case 1 Pag: In this case T i, w) (v, i + aw + (1 a)c, s, w) where e a constant. As result
0 2 VT/ (2t eay= = 0, + ad
0.22 (p(v.) VT. (v= )( O, nu.
In total,
k g = tt ( 7o , s aE2)
Case 2, - In this case T (0, w) = (v - Ureeets i, log u - a, w) so that
0, a VT/ (vt: )ayt 0, 2 g...
0, (u(y.) - VT," - O2 71,
and, thus,
c = N. 0.4 =f 70
Note that 2f is an Omstein-Uhlenbeckk process initialized at 2 and with drift and diffusion matrices
(-o"
14
111
Σ (0)
0
0
112)
02
As a result, we can directly compute 0za2 = ((--s)r, This proves that eq. 19 holds. Eq. 20 then follows directly from 6 and the fact that Ex(y) = k
25

---

From this the results of Section 4.4 follow since the terms oey, vanish whenever s is the spike time of neuron / that is not 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 for each combination of neuron 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 バ to only K² (which is needed anyway to store the weight matrix). In general, the gradient paths can be approximated by simply omitting the terms Ow,y 8.
E Experiments
E.1 Input current estimation
For each combination of sample size and sample 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 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 batch of the same size the sample size and use de 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 out test set of the same size as the sample. We use the RMSProp algorithm with 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 MI 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 stwr://eeenneeeeeoonnomomm refer to the notebook smaxx/notebooks//ttte ssigolloss/sssssoooooeooeeeeeeee neuron. ipynb in the supplementary material.
hold
E.2 Synaptic weight estimation
As above, for each sample size D (256, 512, 1024) we sample data set of spike trains using Algorithm D
with I' = 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:
Ureset = 1.2
A(v) = exp(5(e - 1)
/ (6
a I//4
For each sample size the data was generated using the same randomly sampled weight matrix w which represents a feed-forward network the dimensions described in Section | and which was constructed as follows: for the weight matrix from layer 1 to layer 1 + 1, say w'. we sample each entry from Unif([0.5, 1.5)) and then normalize by 3/ K, where K2 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 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 dj on generated batch and batch sampled from the data set each step. In order to avoid local minimums 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 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/ fnnax/notebooks/ http:///ttbbooossskiinnnnggee wwwissigieeuurreeeeee www.eeeerallaaaa net. ipynb in the supplementary material.
'For updated version of some of the experiments and general implementation of SSNNs in diffrax refer to witggh..coomnnn wwwobb/oo//ocoo//mnmommca
Note that the loss landscape is inherently discontinuous since whenever the parameters are altered in such way that additional spike appears, the expected signature will jump.
26

---

Test loss
Wbros
0.38
Sample size
256
0.36
012
40
1024
0.34
30
0.32
0.30
10
250
500
750
1000
1250
1500
250
500
750
1000
1250
1500
Step
Step
Figure 2: We estimate the synaptic weights across three different sample sizes using the signature kernel MMD truneated depth 3 and stochastic gradient descent with batch size of 128. On the left report the loss on a hold out test set. On the right the mean absolute error between the entries of the currently estimated weight matrix Sxtep and the true weight matrix Wtrue-
NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect paper's contributions and scope?
Answer: [Yes]
Justification: The claims are clearly stated in Section 1.1 which also references where to find the key results.
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 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 long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: There 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 not discussed in the paper.
The authors are encouraged to create 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 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 few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
27

---

The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails 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, worse outcome might be that reviewers discover limitations that aren't It acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
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 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 cross- referenced.
All assumptions should be clearly stated 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 short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that 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) D well as the experimental setups. Additionally, the code is provided in the supplementary material with clear instructions on how to setup the correct environment.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not 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 dataset and/or model, the authors should describe the steps taken make their results reproducible 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 a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to model. In general. releasing code and data often one good way to accomplish this, but reproducibility also be provided via detailed
28

---

instructions for how to replicate the results, access to a hosted model (e.g., g., in the case of 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 submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibilityy the case of closed-source models, may be that access to the model limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions 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
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
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 new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
//nips.
ssioooliccy for more details.
The authors should provide instructions 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. only 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) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimiser, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification: We provide the most important aspects in Section 4 The full details are given in Appendix E
29

---

Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
Justification: We did not prioritise running multiple iterations of the experiments since the main contributions are of a more theoretical nature.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative eтror rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
This information is also given in Appendix E
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
30

---

Answer: [Yes]
Justification:
We have reviewed and made sure that our paper conforms with the code of
ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: The work in the paper of a theoretical nature and there are no clear links to immediate applications with broader societal impacts.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups). privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of MI...
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
31

---

We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
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?
Question:
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 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, paperswi thcode. www.coo//aaaaaaa con/datasete has curated licenses for some datasets. Their licensing guide can help determine the license of 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 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 obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. 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 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 crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
32

---

According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/teview based on the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere 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.
33


Reviews:
[{"review_id": "lmKKMV1qo0", "reviewer": "Reviewer_NV97", "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.", "questions": "What is the computational complexity of SDESolveStep and RootFind?", "limitations": "The authors discussed limitations in Section 5.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper proposes a mathematically rigorous framework to model SSNNs as SDEs with event discontinuities and driven by cadlag rough paths. The paper also identifies sufficient conditions for the existence of gradients, which strictly generalizes previous results only considering ODEs.\n\n2. This paper provides the first gradient-based training of a large class of SSNNs with noise processes in both the spike timing and the network’s dynamics.\n\n3. This paper discusses how the results lead to bio-plausible learning algorithms.", "weaknesses": "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.\n\n2. This paper only conducts experiments on two toy settings. It would be better to discuss how to apply the method to more application conditions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for taking the time to read the paper and provide constructive feedback. We hope that the points below in combination with the author rebuttal answers all the questions that the reviewer may still have.\n\n**W**: 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.\n\n**A**: We agree that we could have spent a little more effort in justifying these assumptions and had also done so in an earlier version, but chose not to here due to space constraints. We propose to make the changes suggested in point 1 of the author rebuttal.\n\n\n**W**: This paper only conducts experiments on two toy settings. It would be better to discuss how to apply the method to more application conditions.\n\n**A**: We agree that the two settings considered in the paper are very limited. The main focus of the paper lies in establishing Theorem 3.1 and explaining its relevance in the context of SSNNs. The algorithm that we present based on the discretize-then-optimize approach should be viewed as a first step in applying these results. Finding ways in which to modify this algorithm to allow for better scaling is left open as an avenue to explore in future work. See also points 2 and 3 of the author rebuttal above for more details and with our suggestions.\n\n**Q**: What is the computational complexity of SDESolveStep and RootFind?\n\n**A**: It is hard to state in generality what the computational complexity of the two operations are since they both depend on the vector field and the method employed. For example, a simple Euler step of the solver requires only a single evaluation of the vector field whereas higher order solves would require more complex operations. We can, however, say something more specific on how the EventSDESolve algorithm scales in the case of a SSNN. See points 2 and 3 of our general rebuttal and the suggestion given therein. The main takeaway is that scaling the discretize-then-optimize algorithm up to a high number of neurons is generally difficult. For these tasks it would appear that the optimize-then-discretize approach is more suitable. Due to time constraints, we have not been able to test out this approach, but as mentioned in the paper, the adjoint equation for the gradients follow readily from Theorem 3.1. We can add a section in the appendix detailing the optimize-then-discretize approach as suggested in point 3 of the author rebuttal.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "DKaGYIWUgV", "reviewer": "Reviewer_sP65", "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.", "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.\n\n2. Can the framework be extended to the case of single spike scenario (for instance in [1,2,3])? \n\n[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).\n\n[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.\n\n[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": "There is no potential negative societal impact.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is generally well-written. Despite its technical nature and numerous technical definitions, it is clear that the authors have made efforts to present the material in an accessible way. \n\n2. The framework in the paper identifies sufficient conditions for the existence of pathwise gradients of solution trajectories and event times with respect to network parameters, and deriving a recursive relation for these gradients.", "weaknesses": "I am not an expert in rough path theory and I don't see any major weaknesses.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for taking the time to read the paper and provide constructive feedback. We hope that the points below in combination with the author rebuttal answers all the questions that the reviewer may still have.\n\n**Q**: 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.\n\n**A**: We note that benchmarking is limited by the fact that, to our knowledge, we are the first to consider stochastic SNNs (in the sense that the ODEs are replaced with SDEs). We agree that more exhaustive empirical evaluations would be very valuable and refer to point 2 and 3 of the author rebuttal above.\n\n\n**Q**: Can the framework be extended to the case of single spike scenario (for instance in [1,2,3])?\n\n**A**: Yes, there is nothing prohibiting us from applying the framework to single spike scenarios. In particular: 1) The main theorem simply states the existence of gradients and how one would compute them when having access to exact solutions of the inter-event SDE. Thus, this result also holds for the first spike times of an output layer of neurons. 2) The signature kernel works for any collection of càdlàg paths. In principle one could just take the first spike time of each output neuron and convert these to spike trains to be fed into the signature kernel MMD.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "57WFiAww5x", "reviewer": "Reviewer_ywki", "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.", "questions": "- Please explain some ambiguous sentences and less known facts:\n  - L211: \"If the diffusion term is absent we recover the gradients in [6]\" from which expression to we recover the gradients in [6]?\n  - Give at least a one sentence description of cadlag (right continuous, left limit). Why is it relevant to the modeling of SNNs?\n  - What is a continuous signature kernel? Can you give an intuition without having to read another technical paper?\n- 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.\n- Can you solve at least one practical SNN training problem?\n- How does the training effort (wall time, memory) compare to existing methods? Is this an impediment to test on a benchmark?", "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.", "rating": 6, "confidence": 2, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "The manuscript offers an attractive solution to the problem of solving SNNs, based on rigorous demonstrations. The derivations and description makes too much use of very specific jargon, making it difficult to evaluate with respect to the existing methods. At this point, if there is indeed practical merit in the proposed approach, it is unlikely to have an impact to the field as is. The lack of testing on standard SNN benchmarks doesn't help this problem. I believe these problems could be fixed by making more parallels with existing work and at least ont benchmark that demonstrates the effectiveness of the proposed method. The yin yang dataset in [55] would be a good candidate.\n\n- Mathematically rigorous demonstration of stochastic SNN gradients is novel and relevant. Strickly speaking, existing methods already demonstrated exact gradients in SNNs in the absence of gradients and when events are known [55], but the extension to stochastic inputs and dynamics is new to my knowledge. Note that in this sense, the statement L260 is already no longer popular belief.\n- Connection to three-factor rules in the feed-forward case extends prior work in the field to stochastic models", "weaknesses": "- 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:\n  \"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\". \n Why adversarial? what is a continuous signature kernel? \n Maybe I'm the wrong reader for this paper, but then this would likely apply to most readers interested in training SNNs.\n\n- 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for taking the time to read the paper and provide constructive feedback. We hope that the points below in combination with the author rebuttal answers all the questions that the reviewer may still have.\n\n**W**: 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.\n\n**A**: We tried our best to merge two very distant communities: Those of computational neuroscience and rough analysis. We believe the value of the contribution is also in this attempt. Due to space limits we cannot go into further details, but we can add the following clarification to the camera-ready version: Training a Neural Stochastic Differential Equation (NSDE) involves minimizing the distance between the path-space distribution generated by the SDE and the empirical distribution supported by observed data sample paths. Various training mechanisms have been proposed in the literature. State-of-the-art performance has been achieved by training NSDEs adversarially as Wasserstein-GANs. However, GAN training is notoriously unstable, often suffering from mode collapse and requiring specialized techniques such as weight clipping and gradient penalty. In [51], the authors introduce a novel class of scoring rules based on signature kernels, a type of characteristic kernel on paths, and use them as the objective for training NSDEs non-adversarially.\n\n\n**W**: 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.\n\n**A**: Please see our points 2 and 3 of the author rebuttal above.\n\n\n**Q**: Please explain some ambiguous sentences and less known facts:\n\n- L211: \"If the diffusion term is absent we recover the gradients in [6]\" from which expression to we recover the gradients in [6]?\n**A**: By this we mean that their eq. (9) is equivalent to our eq. (6). We suggest to add the following in Remark 3.1: \"In particular, eq. (6) for \\(n=1\\) is exactly eq. (9) in [6].\" Note, however, that they only consider the derivative of the first event time and they explicitly model the dependence on \\(t\\) in the event function. We chose to assume that the event function only depends on the state \\(y\\). (This is only for notational simplicity since \\(t\\) can easily be included in the state if need be.)\n\n- Give at least a one sentence description of cadlag (right continuous, left limit). Why is it relevant to the modeling of SNNs?\n**A**: We suggest to add above description as a footnote to the first mention of càdlàg in the main body of text. We note that trajectories of (S)SNNs exhibit jumps due to discontinuous event-triggering and therefore are naturally modelled as càdlàg paths. Spike trains can be viewed as càdlàg step functions counting the number of spikes over time, the membrane potential is right continuous and has left limits (at the discontinuity points the left limit is equal to the firing threshold).\n\n- What is a continuous signature kernel? Can you give an intuition without having to read another technical paper?\n**A**: The signature transform can be viewed as feature map on path space defined as an infinite series of iterated integrals. These iterated integrals play the same role as tensor monomials for vectors on $\\mathbb R^d$, so the signature can be thought as the analogue of a Taylor expansion on path space. The (continuous) signature kernel of two paths is defined as an inner product of their signatures. It serves as a natural measure of (dis)similarity between two curves [*].  We propose to add this explanation to Section 4.1. in the camera-ready version of the paper. \n\n[*] Király, Franz J., and Harald Oberhauser. \"Kernels for sequentially ordered data.\" Journal of Machine Learning Research 20.31 (2019): 1-45.\n\n\n**Q**: Definition 3.1 and 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.\n\n**A**: We understand the confusion that this might cause and would like to point to the paragraph beginning on L219 where we explain the connection between SSNNs and Def. 3.1. The reason for choosing this more general framework is because all of the theory that we develop holds in full generality to what we call Event RDEs (as defined in Section A.4 of the Appendix). This includes many more use cases than SNNs and, in particular, is a generalization (and formalization) of the settings discussed in e.g. [6, 25].\n\n\n**Q**: Can you solve at least one practical SNN training problem?\n\n**A**: Please see the response and suggestions in points 2 and 3 of the author rebuttal.\n\n\n**Q**: How does the training effort (wall time, memory) compare to existing methods? Is this an impediment to test on a benchmark?\n\n**A**: We refer to points 2 and 3 of the author rebuttal.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yC5N2Vx9Od", "reviewer": "Reviewer_2eDR", "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. \n\nThe 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.", "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.\n\n**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? \n\n**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? \n\n[1] Mukhoty, Bhaskar, et al. \"Direct training of snn using local zeroth order method.\" Advances in Neural Information Processing Systems 36 (2024).", "limitations": "Yes, authors do discuss some of the limitations of their work.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper uses stochastic leaky-intergrate-and-fire neurons, which have randomness in the evolution of membrane potential as well as spike timings. It uses the mathematical framework of rough path theory to obtain the gradient of the spike timings with respect to network parameters and formally prove the existence of such a gradient in this framework. The paper is well-written and provides intuitive descriptions along with rigorous mathematical proofs. It brings a novel theoretical insight into spiking neural network training.", "weaknesses": "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for taking the time to read the paper and provide constructive feedback. We hope that the points below in combination with the author rebuttal answers all the questions that the reviewer may still have.\n\n**W**: 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.\n\n**A**: We thank the reviewer for the suggestion and agree that this would be a good addition to the paper. See also point 2 of the author rebuttal.\n\n\n**Q**: **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.\n\n**A**: The event and transition functions determine when the discontinuity happens and how the system jumps at the discontinuity respectively. These functions need to be continuous (as functions of the state) in order for the event times and the solution trajectories to be differentiable (wrt. the initial condition). Note that the discontinuity we are modelling is a discontinuity wrt. time, i.e., the state jumps at certain points in time when an event is triggered. The exact timing and size of these jumps, however, is still allowed to be smooth (as a function of the parameters/initial condition).\n\n\n**Q**: **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?\n\n**A**: The noise from the Brownian motion is implicitly part of eqs. (6) and (7) in two ways: 1) Both expressions depend on the solution trajectories $y^{n-1}_t$ which, in turn, are noisy due to the presence of the Brownian motion. 2) $\\partial y^n_t$ can be shown to satisfy another SDE which in general has a diffusion term as well driven by the same Brownian motion.\n\nIn the deterministic limit we recover the exact gradients of the deterministic SNN with firing exactly at the threshold.\n\nIn the absence of a diffusion term, the dynamics reduce to an ODE in which case exact gradients always exist. In the context of our paper, this is because Assumption 3.4 is satisfied in this case. On the other hand, if a diffusion term is present, we require stochastic firing in order to ensure existence of derivatives via Assumption 3.4.\n\n\n**Q**: **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?\n\n**A**: We believe that this question is answered by points 2 and 3 of the author rebuttal.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We would like to thank all the reviewers for taking the time to go through the paper and providing valuable feedback. We agree with many of the points that have been raised and believe that most, if not all, can be accommodated in a camera-ready version. Apart from minor points and clarifications, we found the main three weaknesses to be shared among all reviews:

1. **Assumptions**: We could have spent more time justifying the assumptions needed to prove existence of solutions of Event SDEs as well as for the main theorem. While we briefly discuss all assumptions in the context of SSNNs, we agree that more could be done. In particular, we agree that the first two Assumptions deserve more explanation and would be happy to add this in the final version. We suggest to add the following after L173: "Assumptions 3.1 and 3.2 simply ensure that an event cannot be triggered immediately upon transitioning. This holds in many settings of interest. For example, for the usual deterministic LIF neuron, $\textnormal{im } \mathcal{T} = 0$ and $\ker \mathcal{E} = 1$ and the length of the refractory period is directly linked to $c$ in Assumption 3.1".  We note that these assumptions are standard in the literature concerning stochastic hybrid processes (see, e.g., our references [31, 32]). We also note that all assumptions are discussed in the context of (S)SNNs after Remark 3.2.

2. **Computational complexity**: More clarity is needed regarding the complexity of the provided algorithm and its alternatives. This is a fair point and we believe it can be answered by adding an additional remark after Remark 3.4 stating the following: "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\ge 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."

3. **Examples**: The primary objective of the paper is 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 above. We propose to add an additional section in the Appendix where we give a derivation for the adjoint equations in a SSNN and provide pseudo-code for an optimize-then-discretize algorithm. As a final comment, we note that the provided Algorithm works for any given event SDE.

Meta-review: 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
Decision: Accept (poster)