High-performance deep spiking neural networks with 0.3 spikes per neuron
Communication by rare, binary spikes is a key factor for the energy
efficiency of biological brains. However, it is harder to train
biologically-inspired spiking neural networks (SNNs) than artificial neural
networks (ANNs). This is puzzling given that theoretical results provide exact
mapping algorithms from ANNs to SNNs with time-to-first-spike (TTFS) coding. In
this paper we analyze in theory and simulation the learning dynamics of
TTFS-networks and identify a specific instance of the vanishing-or-exploding
gradient problem. While two choices of SNN mappings solve this problem at
initialization, only the one with a constant slope of the neuron membrane
potential at threshold guarantees the equivalence of the training trajectory
between SNNs and ANNs with rectified linear units. We demonstrate that training
deep SNN models achieves the exact same performance as that of ANNs, surpassing
previous SNNs on image classification datasets such as MNIST/Fashion-MNIST,
CIFAR10/CIFAR100 and PLACES365. Our SNN accomplishes high-performance
classification with less than 0.3 spikes per neuron, lending itself for an
energy-efficient implementation. We show that fine-tuning SNNs with our robust
gradient descent algorithm enables their optimization for hardware
implementations with low latency and resilience to noise and quantization.
Ana Stanojevic, Stanisław Woźniak, Guillaume Bellec, Giovanni Cherubini, Angeliki Pantazi, Wulfram Gerstner
# Are training trajectories of deep single-spike

###### Abstract

Communication by binary and sparse spikes is a key factor for the energy efficiency of biological brains. However, training deep spiking neural networks (SNNs) with backpropagation is harder than with artificial neural networks (ANNs), which is puzzling given that recent theoretical results provide exact mapping algorithms from ReLU to time-to-first-spike (TTFS) SNNs. Building upon these results, we analyze in theory and in simulation the learning dynamics of TTFS-SNNs. Our analysis highlights that even when an SNN can be mapped exactly to a ReLU network, it cannot always be robustly trained by gradient descent. The reason for that is the emergence of a specific instance of the vanishing-or-exploding gradient problem leading to a bias in the gradient descent trajectory in comparison with the equivalent ANN. After identifying this issue we derive a generic solution for the network initialization and SNN parameterization which guarantees that the SNN can be trained as robustly as its ANN counterpart. Our theoretical findings are illustrated in practice on image classification datasets. Our method achieves the same accuracy as deep ConvNets on CIFAR10 and enables fine-tuning on the much larger PLACES365 dataset without loss of accuracy compared to the ANN. We argue that the combined perspective of conversion and fine-tuning with robust gradient descent in SNN will be decisive to optimize SNNs for hardware implementations needing low latency and resilience to noise and quantization.

## 1 Introduction

Similar to the brain, neurons in spiking neural networks (SNNs) communicate via short pulses called spikes - in striking contrast to artificial neural networks (ANNs) where neurons communicate by the exchange of real-valued signals. While ANNs are the basis of modern artificial intelligence (AI) with impressive achievements [1; 2; 3], their high performance on various tasks comes at the expense of high energy consumption [4; 5; 6]. In general, high energy consumption is a challenge in terms of sustainability or deployment in low-power edge devices [7; 8; 9]. Due to their sparse binary communication scheme, SNNs may offer a potential solution by reducing resource usage in the network [10; 11; 12; 13; 14; 15; 16], but these studies have shown that it is difficult to demonstrate working SNNs which perform at the same level as ANNs.

There exist multiple methods to train the parameters of an SNNs with various advantages and drawbacks. Traditionally, SNNs were trained with plasticity rules observed in biology [17; 18] but it appears more efficient to rely on gradient-descent optimization as done in deep learning (see [19; 20; 21; 22; 23; 24] for theoretical relationships between the plasticity rules and gradient descent). One of the most successful training paradigms for SNNs views the spiking neuron as discrete-time recurrent unit with binary activation and uses a pseudo-derivative or surrogate gradient on the backward pass while keeping the strict threshold function in the forward pass [25; 26; 27; 28; 29].
Other approaches [30, 31, 32] either translate ANN activations into SNN spike counts to train the SNN with the ANN gradients, or use temporal coding with a large number of spikes, both of which can jeopardize the energy efficiency of SNNs.

More recently, and in contrast to spike-count measures in neuroscience [33], it was found in sensory brain areas that neurons also encode information in the exact timing of the first spike, i.e. more salient information leads to earlier spikes [34, 35, 36] which in turn leads to a fast response to the stimuli [37, 38]. While it is possible to train temporally coded spiking neural networks where neurons send _multiple spikes_ to transmit information, we focus in this paper on a _time-to-first-spike_ (TTFS) coding scheme [39] where each neuron fires at most a single spike. The goal of the present study is (1) to analyze theoretically why all attempts at training SNNs with TTFS encoding run into difficulties with deep networks (beyond 4 layers) and (2) to provide a solution to these problems.

Related workThere is a long history of implementations of gradient descent in SNNs with TTFS. In [40] the authors used the Spike Response Model [39] and calculated backpropagation gradients with respect to spike timing and parameters. While the paper states that the learning rule contains an approximation, it turns out to be the exact gradient when the number of spikes is fixed to avoid discontinuities, i.e. spikes do not appear or disappear. This approach was rediscovered recently, extended to other neuron models [14, 41, 42, 43, 44] and applied to small machine learning datasets like MNIST and Fashion-MNIST [14, 45, 42, 44] with architectures of 4 hidden layers or less.

A different line of work avoids training altogether and converts directly an ANN into an SNN. Beyond the classic conversion techniques based on rate coding [46, 47], some studies considered the conversion from ANNs to temporally coded SNNs [48, 49, 50, 51]. While most of them relied on an inexact mapping algorithm or unconventional threshold dynamics, it was recently shown that approximation-free conversion from ANN to TTFS-SNN is possible [51]. However, these results have not shown any benefits outside of the network conversion setting, therefore it remains difficult to understand why gradient descent in the TTFS cannot be used for training or fine-tuning deep SNNs. We build upon this aforementioned work [51] to study the learning dynamics of SNNs.

Contributions of the paperOur work combines the exact backpropagation update steps [41, 42, 43, 44] with an exact mapping between ANNs and SNNs [51], and goes considerably beyond the state-of-the-art with respect to the following points:

* **Even if a TTFS-SNN has an equivalent ReLU ANN, they do not necessarily follow equivalent learning trajectories.** We extend the theory from [51] to provide a reversed mapping from TTFS-SNN to an equivalent ReLU network. We propose the conditions and an adaptive SNN hyperparameter update rule, which are necessary for the equivalence to hold throughout training. Furthermore, the _linearly mappable condition_ from SNN to ANN parameters is identified as the one that guarantees that gradient descent follows the learning trajectories of the equivalent ReLU ANN.
* **Hard instance of the vanishing-gradient problem.** We identify that naively using ANN weight matrix initialization techniques [52, 53, 54] with TTFS-SNN results in a severe instance of the vanishing-gradient problem. We identify the problem analytically and provide a generic recipe to solve it and initialize TTFS-SNN efficiently.
* **Training SNNs to the state-of-the-art accuracy on large datasets.** All previous learning attempts in the TTFS setting were limited to MNIST or Fashion-MNIST datasets and networks of up to 4 layers [14, 45, 55, 43, 42]. We are the first to train TTFS-SNN on CIFAR10 [56], CIFAR100 [56] and PLACES365 [57]. As predicted by the theory, on all datasets, our TTFS-SNN achieves the exact same performance as a ReLU network with the same architecture.
* **Demonstrating the benefits of training under hardware constraints.** Since SNNs can be implemented in low-energy neuromorphic hardware, we test their robustness to quantization of spike times and weights, by fine-tuning the quantized network with our training framework. We show that the latency to decision can be reduced by a factor of four with a performance drop of less than 3 percent on CIFAR10.
## 2 Definition and properties of time-to-first-spike networks

In the following section, we will analyze the gradient descent dynamics in a TTFS setting. To avoid any approximation we follow [51] and study deep spiking neural networks consisting of neurons with triangular post-synaptic integration filters. This model should be viewed as the linearization of a more classical double-exponential filter [14; 45] where all the spikes of a layer arrive within a time window that is small compared to the two time constants of the double-exponential. As we will see in the following section, the theory becomes rigorous under this linearized approximation - we discuss the extension to a spike response model with a double-exponential post-synaptic potential filter in Appendix A.

A time-to-first-spike (TTFS) network modelThe neurons are arranged in \(N\) hidden layers where the spikes of neurons in layer \(n\) are sent to neurons in layer \(n+1\). The layers are either fully-connected (i.e. each neuron receives input from all neurons in the previous layer) or convolutional (i.e. connections are limited to be local and share weights). All connections are feed-forward, i.e. there are no recurrent connections.

_Input layer:_ At the first layer, the analog input to the network represents, for example, the pixel intensity. We assume that the input is scaled to the interval \([0,1]\) and the values are encoded with TTFS coding (a high input pixel intensity \(x_{j}^{(0)}\) leads to an early spike at \(t_{j}^{(0)}\)):

\[t_{j}^{(0)}=\tau_{c}[1-x_{j}^{(0)}]=t_{\max}^{(0)}-\tau_{c}x_{j}^{(0)} \tag{1}\]

where spiking time \(t_{j}^{(0)}\) of neuron \(j\) in the input layer encodes the real-valued input \(x_{j}^{(0)}\in[0,1]\) and \(t_{\max}^{(0)}=\tau_{c}\) is the last possible spike time in the input layer. The conversion parameter \(\tau_{c}\) translates unit-free inputs into time units. In biology \(\tau_{c}\) in sensory areas is in the range of a few milliseconds [34; 35] whereas in hardware devices it could be in the range of microseconds or even shorter.

Neuron dynamicsIn the hidden layers and similarly to [51], the potential \(V_{i}(t)\) of neuron \(i\) models an integrate-and-fire neuron. Given the spike times \(t_{j}^{(n-1)}\) of neurons \(j\) in the previous layer, the potential \(V_{i}^{(n)}\) of neuron \(i\) in layer \(n\) follows the dynamics:

\[\tau_{c}\frac{\mathrm{d}V_{i}^{(n)}}{\mathrm{d}t}=\alpha_{i}^{(n)}H(t-t_{\text {min}}^{(n)})+\sum_{j}W_{ij}^{(n)}H(t-t_{j}^{(n-1)}) \tag{2}\]

where \(t_{\min}^{(n)}\) is a constant which will become by construction a lower bound of the earliest possible spike time in layer \(n\), \(\alpha_{i}^{(n)}\) is a positive scalar which can be seen as the weight of an external spike at \(t_{\min}^{(n)}\), \(W_{ij}^{(n)}\) is the synapse strength from neuron \(j\) to neuron \(i\), and \(H\) denotes the Heaviside function

Figure 1: **Network of TTFS neurons.****a.** A feed-forward TTFS architecture. **b.** Potential \(V_{i}^{(n)}\) for different neurons \(i\) as a function of time. The intial slope is always zero. For \(\alpha_{i}^{(n)}=1\) slopes of all neurons \(i\) are increased by \(1\) at \(t_{\text{max}}^{(n-1)}=t_{\text{min}}^{(n)}\). **c.** For \(\alpha_{i}^{(n)}=1-\sum_{j}W_{ij}^{(n)}\) slopes of all neurons \(i\) are equal to \(1\) for \(t\geq t_{\text{min}}^{(n)}\).
which takes a value of \(1\) for positive inputs and is \(0\) otherwise. When the potential \(V_{i}^{(n)}\) reaches the threshold \(\vartheta_{i}^{(n)}\), neuron \(i\) generates a spike at time \(t_{i}^{(n)}\) and sends it to the next layer. Once a neuron spikes we assume a very large refractory period to ensure that every neuron spikes at most once. According to Eq. (2), each input spike \(t_{j}^{(n-1)}\) changes the _slope_ of the potential by a fixed amount proportional to \(W_{ij}^{(n)}\). Integration of Eq. (2) therefore leads to a piecewise linear behaviour of the potential, see Fig. 1. Without loss of generality, we assume \(V\) to be unit-free and so are the parameters \(W_{ij}^{(n)}\) and \(\alpha_{i}^{(n)}\), whereas \(t\) and \(\tau_{c}\) have units of time. Rescaling time by \(t\rightarrow(t/\tau_{c})\) would remove the units, but we keep it in the equations to show the role of the conversion factor \(\tau_{c}\). Note that this is the same spiking neuron model as defined in [51] with the minor modification that the ramping input of strength \(\alpha_{i}^{(n)}\) arrives at \(t_{\text{min}}^{(n)}\) and not \(t_{\text{min}}^{(n-1)}\), which will simplify our equations in the following. We initialize \(\vartheta_{i}^{(n)},t_{\text{min}}^{(n)},t_{\text{max}}^{(n)}\) so that all the neurons of layer \(n\) spike once in the interval \([t_{\text{min}}^{(n)},t_{\text{max}}^{(n)}]\). The threshold \(\vartheta_{i}^{(n)}\) is defined as \(\vartheta_{i}^{(n)}\stackrel{{\text{def}}}{{=}}\tilde{\vartheta }_{i}^{(n)}+D_{i}^{(n)}\) where \(D_{i}^{(n)}\) is a model parameter initialized at 0 and \(\tilde{\vartheta}_{i}^{(n)}\) is the base threshold. We define a maximum spike time \(t_{\text{max}}^{(n)}\), after which emission of a spike is forced in all neurons of layer \(n\) which have not spiked yet. The construction of \(t_{\text{min}}^{(n)}\) and \(t_{\text{max}}^{(n)}\) is recursive as \(t_{\text{min}}^{(n)}\stackrel{{\text{def}}}{{=}}t_{\text{max}}^ {(n-1)}\). Details of the choice of the base threshold and \(t_{\text{max}}^{(n)}\) are given Appendix B.

Adaptive \(t_{\text{max}}^{(n)}\) parametersDuring training as we update the network parameters \(W_{ij}^{(n)}\) and \(D_{i}^{(n)}\) the hyperparameters like \(t_{\text{max}}^{(n)}\) need to be changed as well such that the condition that all the neurons of layer \(n\) spike once in the interval \([t_{\text{min}}^{(n)},t_{\text{max}}^{(n)}]\) remains true. We suggest a new adaptive update rule which moves \(t_{\text{max}}^{(n)}\). Note that this is an addition to the model in [51] where an adaptive \(t_{\text{max}}^{(n)}\) was not necessary since the parameters were fixed. Formally, when processing the training dataset, we update \(t_{\text{max}}^{(n)}\) as follows:

\[\Delta t_{\text{max}}^{(n)}=\begin{cases}\gamma(t_{\text{max}}^{(n)}-\text{ min}_{i,\mu}t_{i}^{(n)})-(t_{\text{max}}^{(n)}-t_{\text{min}}^{(n)}),&\text{if }t_{\text{max}}^{(n)}-t_{\text{min}}^{(n)}<\gamma(t_{\text{max}}^{(n)}-\text{ min}_{i,\mu}t_{i}^{(n)})\\ 0,&\text{otherwise}\end{cases} \tag{3}\]

The minimum operator iterates over all neurons \(i\) and input samples \(\mu\) in the batch and \(\gamma\) is a constant. After this update, we change the subsequent time window accordingly so that: \(t_{\text{min}}^{(n+1)}=t_{\text{max}}^{(n)}\), and we iterate over all layers sequentially. The base threshold \(\tilde{\vartheta}_{i}^{(n)}\) is then updated accordingly, see Appendix B. This adaptation effectively moves all the spikes \(t_{j}^{(n)}\) away from the boundary \(t_{\text{min}}^{(n)}\). For simplicity, in the theory section, we consider that this update has reached an equilibrium, so we consider that \(t_{\text{max}}^{(n)}\), \(t_{\text{min}}^{(n)}\), and \(\tilde{\vartheta}_{i}^{(n)}\) are constants w.r.t. the SNN parameters, the condition \(t_{\text{min}}^{(n+1)}=t_{\text{max}}^{(n)}\) is always satisfied and all the spikes of layer \(n\) arrive within \([t_{\text{min}}^{(n)},t_{\text{max}}^{(n)}]\).

Output layerThe output layer has index \(N+1\) and contains non-spiking read-out neurons. Each neuron \(m\) simply integrates input spikes coming from layer \(N\) without firing. Integration of \(V_{m}^{(N+1)}\) stops at time \(t_{\text{min}}^{(N+1)}\) and the _softmax_ and the standard cross-entropy loss \(\mathcal{L}\) are calculated using real-valued potentials, analogous to real-valued activations of ANNs.

General reverse mapping from SNN to ANNBuilding upon the conversion method from ANN with rectified linear units (ReLUs) to TTFS networks [51], we now describe a reversed mapping strategy defining uniquely the parameters of an equivalent ANN for given SNN parameters. This mapping will be a fundamental pillar in the following theoretical analysis. In the most general case, to find a ReLU network with weights \(w\) and bias \(\mathbf{b}\) which is equivalent to our SNN model, we define:

\[B_{i}^{(n)}\stackrel{{\text{def}}}{{=}}\alpha_{i}^{(n)}+\sum_{k}W_{ ik}^{(n)}\quad\text{and}\quad\ w_{ij}^{(n)}\stackrel{{\text{def}}}{{=}} \frac{W_{ij}^{(n)}}{B_{i}^{(n)}}\quad\text{and}\quad\ b_{i}^{(n)}\stackrel{{ \text{def}}}{{=}}-\frac{\vartheta_{i}^{(n)}}{B_{i}^{(n)}}+\frac{t_{\text{ max}}^{(n)}-t_{\text{min}}^{(n)}}{\tau_{c}}. \tag{4}\]

Here, \(B_{i}^{(n)}\) has a simple interpretation: it is the slope of the potential at the moment of threshold crossing in the SNN if time is measured in units of \(\tau_{c}\) (see Eq. (2)). We call \(B_{i}^{(n)}\) the'slope-at-threshold factor' and it will play an important role in the following. Then if we define the ANN
activation at the input layer as the pixel intensity \(\mathbf{x}^{(0)}\), Eq. (4) defines uniquely a ReLU network with activations \(\mathbf{x}^{(n)}\) such that (this is the reciprocal mapping inspired by [51]):

\[\mathbf{x}^{(n)}\tau_{c}=t_{\text{max}}^{(n)}-\mathbf{t}^{(n)}. \tag{5}\]

At the output layer, we resort directly to the simpler parameter mapping of the output layer from [51]: with \(w_{ij}^{(N+1)}\stackrel{{\text{\tiny def}}}{{=}}W_{ij}^{(N+1)}\) and \(b_{i}^{(N+1)}\stackrel{{\text{\tiny def}}}{{=}}\alpha_{i}^{(N+1)} (t_{\text{max}}^{(N)}-t_{\text{min}}^{(n)})\) the logits and the cross-entropy loss \(\mathcal{L}\) are also equal in the SNN and the equivalent ANN.

Proof.: Starting from the introduced SNN definition, we compute analytically the spikes at time \(t_{i}^{(n)}\) in the SNN. In case the potential \(V_{i}^{(n)}\) reaches the threshold \(\vartheta_{i}^{(n)}\) before \(t_{\text{max}}^{(n)}\), the spiking condition \(\vartheta_{i}^{(n)}=V_{i}^{(n)}(t_{i}^{(n)})\) yields:

\[\tau_{c}\vartheta_{i}^{(n)}=\alpha_{i}^{(n)}(t_{i}^{(n)}-t_{\text{min}}^{(n)} )H(t_{i}^{(n)}-t_{\text{min}}^{(n)})+\sum_{j}W_{ij}^{(n)}(t_{i}^{(n)}-t_{j}^{(n -1)})H(t_{i}^{(n)}-t_{j}^{(n-1)}). \tag{6}\]

Since we constructed the SNN such that \(t_{i}^{(n)}\) arrives in the time window \([t_{\text{min}}^{(n)},t_{\text{max}}^{(n)}]\), all the terms \(H(\cdot)\) are equal to \(1\) in the previous equation. So for any spike \(t_{i}^{(n)}\) arriving before \(t_{\text{max}}^{(n)}\), we have:

\[t_{i}^{(n)}=\frac{\tau_{c}\vartheta_{i}^{(n)}+\alpha_{i}^{(n)}t_{\text{min}}^{ (n)}+\sum_{j}W_{ij}^{(n)}t_{j}^{(n-1)}}{\alpha_{i}^{(n)}+\sum_{k}W_{ik}^{(n)}} =\frac{A_{i}^{(n)}}{B_{i}^{(n)}}\;. \tag{7}\]

We can already identify the slope-at-threshold \(B_{i}\) in the denominator, then we replace \(\alpha_{i}^{(n)}\) by \(B_{i}^{(n)}-\sum_{k}W_{ij}^{(n)}\) in \(A_{i}^{(n)}\) and subtract \(t_{\text{max}}^{(n)}\) on both sides:

\[t_{i}^{(n)}-t_{\text{max}}^{(n)}=\tau_{c}\frac{\vartheta_{i}^{(n)}}{B_{i}^{(n )}}+t_{\text{min}}^{(n)}-t_{\text{max}}^{(n)}+\sum_{j}\frac{W_{ij}^{(n)}}{B_{i }^{(n)}}(t_{j}^{(n-1)}-t_{\text{min}}^{(n)}). \tag{8}\]

Using this identity, and using that the rectified linear unit is in its operating regime \(x_{i}>0\) if and only if the spiking neuron \(i\) fires before \(t_{\text{max}}\), one can now prove by induction that the definition Eq. (4) defines an equivalent ReLU network satisfying the identity from Eq. (5). 

## 3 Analysis of learning dynamics

The linearly mappable conditionWe now define a specific choice of SNN, so-called _linearly mappable SNN_ and we will show that it satisfies the theoretical conditions for robust SNN training via gradient descent optimization. _The linearly mappable condition_ is defined by the choice of \(\alpha_{i}^{(n)}\):

\[\alpha_{i}^{(n)}=1-\sum_{j}W_{ij}^{(n)}. \tag{9}\]

This choice implies that the slope-at-threshold \(B_{i}^{(n)}=1\), and results in the linear mapping formula:

\[w_{ij}^{(n)}\stackrel{{\text{\tiny def}}}{{=}}W_{ij}^{(n)}\quad \text{ and }\quad b_{i}^{(n)}\stackrel{{\text{\tiny def}}}{{=}}- \vartheta_{i}^{(n)}+\frac{t_{\text{max}}^{(n)}-t_{\text{min}}^{(n)}}{\tau_{c}} \tag{10}\]

Vanishing-gradient problem for deep SNNsPrevious TTFS-SNNs with exact gradients use mostly shallow networks containing one hidden layer [12; 45; 43; 44], or at most 4 hidden layers [42]. The question arises why the exact gradient approach does not scale to larger networks. We demonstrate in this section that TTFS networks are prone to yield vanishing or exploding gradients (vanishing-gradient problem [54; 52]). To solve this problem, we show that a tight balance has to be respected between the weight initialization and the slope-at-threshold vector \(\textbf{B}^{(n)}\). This analysis will result in the definition of robust initialization schemes for TTFS networks.

The vanishing-gradient problem has been studied exhaustively in ANNs [52]. Similarly, TTFS networks are also subject to this problem. To see this, one has to observe that the network state at
layer \(n\) is summarized by the vector of spike timings \(\mathbf{t}^{(n)}\) such that the loss with respect to the weight parameter at layer \(n\) factorizes as:

\[\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\mathbf{W}^{(n)}}=\frac{\mathrm{d} \mathcal{L}}{\mathrm{d}\mathbf{V}^{(N+1)}}\frac{\mathrm{d}\mathbf{V}^{(N+1)}}{ \mathrm{d}\mathbf{t}^{(N)}}\frac{\mathrm{d}\mathbf{t}^{(N)}}{\mathrm{d} \mathbf{t}^{(N-1)}}\cdots\frac{\mathrm{d}\mathbf{t}^{(n+1)}}{\mathrm{d} \mathbf{t}^{(n)}}\frac{\mathrm{d}\mathbf{t}^{(n)}}{\mathrm{d}\mathbf{W}^{(n)}}. \tag{11}\]

where \(\mathbf{V}^{(N+1)}\) is a vector containing potentials of neurons in layer \(N+1\) at time \(t_{\text{min}}^{(N+1)}\). Hence, if the product of Jacobians: \(\frac{\mathrm{d}\mathbf{t}^{(n+1)}}{\mathrm{d}\mathbf{t}^{(n)}}\) is naively defined, the amplitude of this gradient might vanish or explode exponentially fast as the number of layers becomes large. As analyzed in [52, 58], a way to solve this problem is to make sure that the largest eigenvalues of the Jacobian \(\frac{\mathrm{d}\mathbf{t}^{(n+1)}}{\mathrm{d}\mathbf{t}^{(n)}}\) are close to \(1\) in absolute value.

We now compute analytically the Jacobian of the SNN. It requires the definition \(M_{i}^{(n)}\) which is \(1\) if and only if spike \(t_{i}^{(n)}\) arrives before \(t_{\max}^{(n)}\). We also denote with \(M^{(n)}\) the matrix containing \(M_{i}^{(n)}\) on the diagonal and \(0\) elsewhere, with \(\mathbf{t}^{(n)}\) a vector of spike times in layer \(n\) and with \(\mathrm{B}^{(n)}\) a diagonal matrix of slope-at-threshold factors. By deriving Eq. (7) we find that the Jacobian of the network can be written as (\(\cdot\) is the matrix multiplication):

\[\frac{\mathrm{d}\mathbf{t}^{(n)}}{\mathrm{d}\mathbf{t}^{(n-1)}}=M^{(n)}\cdot \frac{1}{B^{(n)}}\cdot W^{(n)} \tag{12}\]

We can now analyze the conditions for which the vanishing-or-exploding gradient problems are solved. We observe primarily that (1) the eigenvalues of this Jacobian are strongly determined by the slope-at-threshold \(B^{(n)}\) and not only by the weight matrix \(W^{(n)}\) as in ANN; then (2) the eigenvalues of the Jacobian of the SNN are the same as the eigenvalues of the Jacobian of the equivalent ANN. To see this, one may recall from Eq. (4) that the ANN weights are \(w^{(n)}=\frac{1}{B^{(n)}}\cdot W^{(n)}\) and \(M_{i}^{(n)}\) is 1 if and only if the equivalent ReLU unit is saturated.

Robust initialization of TTFS networksBefore defining a generic recipe for initializing the weight matrix \(W^{(n)}\), we illustrate why using naively the standard deep learning recipes with SNN results in vanishing or exploding gradients. In Fig. 2, we demonstrate numerically that this naive approach faces the vanishing-gradient problem. We initialized the weight matrix of an SNN with \(W^{(n)}=\frac{1}{\sqrt{340}}\mathcal{N}(0,1)\) where \(340\) is the number of units in the layers (this is one of the many standard choices in deep learning) so the eigenvalue of \(W^{(n)}\) with largest absolute value is close to \(1\). We can use this matrix to estimate the eigenvalue spectrum of the SNN at initialization. Following classical work in the ANN literature [58, 59], we assume that \(M^{(n)}\) has a small impact on the distribution of the eigenvalues, and we can display the eigenvalue spectrum of \(w^{(n)}=\frac{1}{B^{(n)}}\cdot W^{(n)}\) which are closely related to the eigenvalue spectrum of the SNN Jacobian \(\frac{\mathrm{d}\mathbf{t}^{(n)}}{\mathrm{d}\mathbf{t}^{(n-1)}}\). As shown in Fig. 2a, this naive initialization produces multiple eigenvalues with modulus larger than \(1\) when \(\alpha_{i}^{(n)}=1\); this eigenvalue spectrum leads to an explosion of the gradient norm in backpropagation.

Our theory also provides a recipe to initialize the SNN outside of the _linearly mappable condition_. Since we know from Eq. (12) that a good SNN initialization requires the equivalent ANN to have a

Figure 2: **Eigenvalues of the SNN Jacobian under naive initialization and solution.****a.** In the general case, when the SNN is initialized with standard deep learning initialization, the eigenvalues of the SNN spread beyond the unit circle. **b.** It is possible to correct for this using the inverse relation between \(w^{(n)}\) and \(W^{(n)}\). With the linearly mappable condition, this correction is not necessary.
good initialization, we can first choose the matrix in the ANN parameter space, and map it to the SNN initialization with the inverse map of Eq. (4). This is for instance automatically corrected with the _linearly mappable condition_ in Fig. 2b since the SNN and ANN are the same, so the eigenvalues stay tightly within the unit circle, showing numerically that the vanishing-gradient problem is avoided.

Biased gradient descent trajectory with the generic mappingBeyond initialization, we also analyze whether the gradient descent trajectory in the SNN parameter space necessarily follows the gradient descent trajectory of the equivalent ReLU network. To describe the gradient descent trajectory of the SNN, we consider a gradient descent step with learning rate \(\eta\) when applying backpropagation to the SNN: \(\Delta W_{ij}^{(n)}=-\eta\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}W_{ij}^{(n)}}\), and compute the corresponding update in the space of the ANN parameters. We denote with \(\delta w_{ij}^{(n)}\) the update in ANN parameter space, and we use \(\mathcal{M}_{\mathbf{\alpha}}\) to denote the mapping formula such that \(w_{ij}^{(n)}=\mathcal{M}_{\mathbf{\alpha}}(W_{ij}^{(n)})\); see Eq. (4). If we assume that \(\alpha_{i}^{(n)}\) is a constant independent of \(W_{ij}^{(n)}\), we find: \(\delta w_{ij}^{(n)}=\mathcal{M}_{\mathbf{\alpha}}(W_{ij}^{(n)}-\eta\frac{d \mathcal{L}}{dW_{ij}^{(n)}})-\mathcal{M}_{\mathbf{\alpha}}(W_{ij}^{(n)})\). Assuming a small learning rate, we make a first-order approximation using the derivative \(\frac{\mathrm{d}\mathcal{M}_{\mathbf{\alpha}}}{\mathrm{d}W_{ij}^{(n)}}=\frac{ \mathrm{d}w_{ij}^{(n)}}{\mathrm{d}W_{ij}^{(n)}}=\frac{B_{ij}^{(n)}-W_{ij}^{(n )}}{(B_{i}^{(n)})^{2}}\) leading to the approximate update in ANN parameter space:

\[\delta w_{ij}^{(n)}\approx-\eta\,\frac{\mathrm{d}\mathcal{M}_{\mathbf{\alpha}}}{ \mathrm{d}W_{ij}^{(n)}}\,\delta W_{ij}^{(n)}=-\eta\frac{\mathrm{d}\mathcal{M} _{\mathbf{\alpha}}}{\mathrm{d}W_{ij}^{(n)}}\,\frac{\mathrm{d}\mathcal{L}}{ \mathrm{d}W_{ij}^{(n)}}=-\eta\left[\frac{\mathrm{d}w_{ij}^{(n)}}{\mathrm{d}W _{ij}^{(n)}}\right]^{2}\,\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}w_{ij}^{(n)}} \tag{13}\]

The difference between Eq. (13) and a direct ANN update obtained through gradient descent \(\delta w_{ij}^{(n)}\propto\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}w_{ij}^{(n)}}\) cannot be corrected with a different learning rate \(\eta\), because the multiplicative bias in Eq. (13) changes for every neuron pair \((i,j)\) and algorithmic iteration. We conclude that in general, the gradient descent trajectory in SNN is "biased", meaning that it is impossible to find a naive gradient descent trajectory in the ANN for which SNN and ANN remain equivalent from initialization to convergence. Under the linearly mappable condition, the multiplicative bias disappears since \(\frac{\mathrm{d}w_{ij}^{(n)}}{\mathrm{d}W_{ij}^{(n)}}=1\) and this is the choice made in Section 4. An alternative might be to work with specifically designed'metrics' [61] that counterbalances the multiplicative factor in Eq. (13).

The difficulty to train an SNN with \(\alpha_{i}^{(n)}=1\) is illustrated in Fig. 3. Both SNNs (with or without the linearly mappable condition) are initialized to be equivalent to the same ReLU ANN which solves the vanishing-gradient problem at initialization. Nevertheless, only the SNN with the linearly mappable condition follows the ANN whereas the other one diverges away from the true ReLU trajectories after \(20\) epochs (Fig. 3b). This is true despite using the small learning rate.

Figure 3: **Learning trajectories are biased without the linearly mappable condition.****a.** Training of an 8-layer SNN on the MNIST dataset [60] with and without the linearly mappable condition. Both are initialized correctly: we make sure the eigenvalues of the Jacobian at initialization lie within the unit circle. The SNN (light blue) follows the same training curve as the equivalent ANN (red-dashed) only under the linearly mappable condition. **b.** Weights \(W^{(n)}\) of the SNN and \(w^{(n)}\) of the equivalent ANN remain the same (Cosine Similarity=1) throughout training under the linearly mappable condition, otherwise \(\frac{W^{(n)}}{\mathrm{B}^{(n)}}\) is the same as \(w^{(n)}\) in the beginning but diverges during training.
## 4 Benchmark results

In the following we always consider SNN initialized with the linearly mappable condition and trained with the Adam optimizer and exponential learning rate schedule (See Appendix C simulation details).

Comparison with previous TTFS trainingWe report the performance on the MNIST [60] and Fashion-MNIST [62] to compare with previous implementation of TTFS training paradigms. We tested a 16-layer fully-connected SNN and a ConvNet SNN (similar to LeNet5). As expected from the theory our network achieves the same performance as the ReLU ANN as seen in Figure 4 and Table 1. The performance is therefore better than all previous TTFS implementations which were limited to \(4\) layers (see Appendix C).

State-of-the art TTFS performance on CIFAR100 and PLACES365Previously, tackling larger-scale image datasets like CIFAR100 [56] or PLACES365 [57] (similarly large as ImageNet) was considered impossible. We propose to combine conversion from pre-trained VGG16 and fine-tuning with our approach to build competitive TTFS ConvNet. In Fig. 4b we use a pre-trained VGG16 architecture downloaded from an online repository [63; 64] and map it to the SNN without loss of performance (similarly to [51]). The networks are then fine-tuned with a reduced learning rate. In Table 2 we see results for different datasets, fine-tuning provides an increase of the SNN accuracy by \(1.76\%\) on CIFAR100 and \(1.17\%\) on PLACES365. We are not aware of any TTFS method achieving higher accuracies. More interestingly, fine-tuning SNNs promises to be most useful when the SNN performance is degraded through conversion, for instance, because of hardware constrains like quantization as demonstrated in the next section.

\begin{table}
\begin{tabular}{l c c} Dataset & \multicolumn{2}{c}{Acc [\%]} \\  & ReLU & SNN \\ \hline \hline MNIST & 99.57 \(\pm\) 0.01 & 99.57 \(\pm\) 0.00 \\ f-MNIST & 94.24 \(\pm\) 0.02 & 94.26 \(\pm\) 0.03 \\ CIFAR10 & 93.68 \(\pm\) 0.02 & 93.69 \(\pm\) 0.001 \\ \hline \end{tabular}
\end{table}
Table 1: Resulting when training TTFS-SNNs on MNIST, f-MNIST and CIFAR10

\begin{table}
\begin{tabular}{l c c c c c c} Dataset & Image size & Classes & Acc [\%] & w/o FT & Acc [\%] & w/ FT [\%] \\  & & & ReLU & SNN & ReLU & SNN \\ \hline \hline CIFAR100 & 32 \(\times\) 32 \(\times\) 3 & 100 & 70.48 & 70.48 & 72.23 \(\pm\) 0.06 & 72.24 \(\pm\) 0.06 \\ PLACES365 & 224 \(\times\) 224 \(\times\) 1 & 365 & 52.69 & 52.69 & 53.86 \(\pm\) 0.01 & 53.86 \(\pm\) 0.02 \\ \hline \end{tabular}
\end{table}
Table 2: VGG16 architecture before and after fine-tuning, for our SNN and the equivalent ANN

Figure 4: **Performance of the SNN under the linearly mappable condition a.** Training deep SNN and ReLU network yields same average performance. **b.** Fine-tuning (FT) of VGG16 SNN network which is initialized with the weights of a pre-trained VGG16 ReLU network. Table 2 shows results after fine-tuning for both networks.
Mitigating quantization, noise effects and reducing latencyLet's consider a scenario in which a ReLU network was pre-trained with full precision weights. After mapping to the SNN it is assumed to be deployed on a device with parameter noise, limited temporal resolution or limited weight precision. We use fine-tuning of the VGG16 to recover the SNN accuracy in all these situations in Fig. 5. All the experiments are done on CIFAR10 dataset with 10 epochs of fine-tuning. In all three cases (spike time jitter, time-step quantization or SNN weight quantization) fine-tuning enables to recover the performance. We demonstrate in particular TTFS VGG16 networks achieving higher than \(90\%\) accuracy on CIFAR 10 with 16 time-steps per layer or weights quantized on \(4\) bits. We also investigated whether it is possible to improve the classification latency through fine-tuning by reducing the intervals \([t_{\text{min}}^{(n)},t_{\text{max}}^{(n)}]\) after conversion from ANN. Doing this we indeed improve latency but the SNN diverges away from the pre-trained ANN. Through fine-tuning, performance higher than \(90\%\) test accuracy is recovered, even when the latency is improved by a factor of \(4\).

## 5 Discussion and future work

In this work we solved the hard instance of a vanishing-gradient problem for single-spike neural networks. Moreover, we showed that through application of the linear-mapping condition the learning trajectories of ReLU ANN and SNN become equivalent. Based on this result, we demonstrated, to the best of our knowledge for the first time, that training of deep single-spike neural networks with sixteen layers yields identical performance as ReLU ANN on large datasets such CIFAR100 and PLACES365. In the future we plan to train even deeper networks and more sophisticated architectures such as ResNets, but this requires to map skip connections to the SNN model which is not trivial.

The work will be probably most impactful when implemented in either digital or analog SNN hardware. We are able to fine-tune the single-spike neural network to adapt the SNN to specific device constraints. Moreover, the learning can be generalized to leaky neuronal dynamics which further addresses the imperfections in hardware elements. After downloading the pre-trained network on the device we envision a continual online learning on chip with energy efficient and low latency inference. Our demonstration that gradient descent is possible in deep SNNs might provide in the future the opportunities to derive a completely local, hardware-friendly, training algorithm.

Figure 5: **Fine-tuning (FT) VGG16 for quantization/noise robustness or low latency.** In all cases after only 10 epochs of fine-tuning on CIFAR10 (purple) the performance of the initially mapped network (blue) is significantly improved. **a.** Adding noise with certain standard deviation (SD) to all spiking times in the network. **b.** Quantizing spiking times in the network to a given number of time steps per layer. **c.** Representing all weights \(W_{ij}^{(n)}\) with given number of bits. **d.** Reducing the latency by reducing the size of \([t_{\text{min}}^{(n)},t_{\text{max}}^{(n)}]\).